{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.支持向量机（SVM）和SMO算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import *\n",
    "from time import sleep\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1.简化版SMO算法处理小规模数据集\n",
    "### 6.1.1.辅助函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def loadDataSet(fileName):\n",
    "    \"\"\"\n",
    "    读取数据集\n",
    "    参数：\n",
    "        fileName -- 文件名\n",
    "    返回：\n",
    "        dataMat -- 数据矩阵\n",
    "        label -- 标签矩阵\n",
    "    \"\"\"\n",
    "    #新建数据和标签矩阵列表\n",
    "    dataMat = []; labelMat = []\n",
    "    #获得文件指针\n",
    "    fr = open(fileName)\n",
    "    #一行行读取\n",
    "    for line in fr.readlines():\n",
    "        #以制表符分割\n",
    "        lineArr = line.strip().split('\\t')\n",
    "        #数据1、2列添加到数据矩阵\n",
    "        dataMat.append([np.float(lineArr[0]), np.float(lineArr[1])])\n",
    "        #数据3列添加到标签矩阵\n",
    "        labelMat.append(np.float(lineArr[2]))\n",
    "    return dataMat,labelMat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def selectJrand(i,m):\n",
    "    \"\"\"\n",
    "    随机选取和i不同的j\n",
    "    参数:\n",
    "        i -- 第一个alpha的下标\n",
    "        m -- alpha的个数\n",
    "    返回：\n",
    "        j -- 随机的值\n",
    "    \"\"\"\n",
    "    j=i\n",
    "    while (j==i):\n",
    "        j = np.int(np.random.uniform(0,m))\n",
    "    return j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def clipAlpha(aj,H,L):\n",
    "    \"\"\"\n",
    "    约束alpha\n",
    "    参数：\n",
    "        aj -- 迭代得到的未约束alpha的值\n",
    "        H -- 上界\n",
    "        L -- 下界\n",
    "    返回：\n",
    "        aj -- 约束后的alpha\n",
    "    \"\"\"\n",
    "    if aj > H: \n",
    "        aj = H\n",
    "    if L > aj:\n",
    "        aj = L\n",
    "    return aj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1.2.简化的SMO算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def smoSimple(dataMatIn, classLabels, C, toler, maxIter):\n",
    "    \"\"\"\n",
    "    简化的SMO算法\n",
    "    参数：\n",
    "        dataMatIn -- 输入数据矩阵\n",
    "        classLabels -- 标签类列表\n",
    "        C -- 权重超参数\n",
    "        toler -- 容错率\n",
    "        maxIter -- 最大迭代次数\n",
    "    返回：\n",
    "        b -- 更新的偏差参数\n",
    "        alphas -- 更新的alpha参数向量\n",
    "    \"\"\"\n",
    "    #数据和标签转化为矩阵形式\n",
    "    dataMatrix = np.array(dataMatIn); labelMat = np.array(classLabels).reshape(-1, 1)\n",
    "    #偏差参数初始化\n",
    "    #记录数据矩阵的维度\n",
    "    b = 0; m,n = np.shape(dataMatrix)\n",
    "    #初始化alpha向量为全0\n",
    "    alphas = np.zeros((m,1))\n",
    "    #初始化迭代次数计数器\n",
    "    iter = 0\n",
    "    #当迭代次数小于最大迭代次数时\n",
    "    while (iter < maxIter):\n",
    "        #alpha对修改次数计数器\n",
    "        alphaPairsChanged = 0\n",
    "        #寻找可优化的alpha\n",
    "        for i in range(m):\n",
    "            fXi = np.float(np.dot((alphas*labelMat).T, np.dot(dataMatrix, dataMatrix[i:i+1,:].T))) + b\n",
    "            Ei = fXi - np.float(labelMat[i, 0])\n",
    "            #检查是否有元素违背KKT条件\n",
    "            if ((labelMat[i, 0]*Ei < -toler) and (alphas[i, 0] < C)) or ((labelMat[i, 0]*Ei > toler) and (alphas[i, 0] > 0)):\n",
    "                #随机选择j\n",
    "                j = selectJrand(i,m)\n",
    "                fXj = np.float(np.dot((alphas*labelMat).T, np.dot(dataMatrix, dataMatrix[j:j+1,:].T))) + b\n",
    "                Ej = fXj - np.float(labelMat[j, 0])\n",
    "                #提取旧的alphai和alphaj\n",
    "                alphaIold = alphas[i, 0].copy(); alphaJold = alphas[j, 0].copy();\n",
    "                #如果i和j的标签不同\n",
    "                if (labelMat[i, 0] != labelMat[j, 0]):\n",
    "                    L = max(0, alphas[j, 0] - alphas[i, 0])\n",
    "                    H = min(C, C + alphas[j, 0] - alphas[i, 0])\n",
    "                else:\n",
    "                    L = max(0, alphas[j, 0] + alphas[i, 0] - C)\n",
    "                    H = min(C, alphas[j, 0] + alphas[i, 0])\n",
    "                #若上界等于下界，就执行下一次循序\n",
    "                if L==H:  continue\n",
    "                #计算最优修改量\n",
    "                eta = 2.0 * np.dot(dataMatrix[i:i+1,:], dataMatrix[j:j+1,:].T) - np.dot(dataMatrix[i:i+1,:], dataMatrix[i:i+1,:].T) - np.dot(dataMatrix[j:j+1,:], dataMatrix[j:j+1,:].T)\n",
    "                if eta >= 0:  continue\n",
    "                #更新alphaj\n",
    "                alphas[j, 0] -= labelMat[j, 0]*(Ei - Ej)/eta\n",
    "                #约束alphaj\n",
    "                alphas[j, 0] = clipAlpha(alphas[j, 0],H,L)\n",
    "                #判断精度是否满足要求，是，则执行下一次循环\n",
    "                if (abs(alphas[j, 0] - alphaJold) < 0.00001):  continue\n",
    "                #修改i的量和j相同，但是方向相反，保证约束条件成立\n",
    "                alphas[i, 0] += labelMat[j, 0]*labelMat[i, 0]*(alphaJold - alphas[j, 0])\n",
    "                #计算新的b\n",
    "                b1 = b - Ei- labelMat[i, 0]*(alphas[i, 0]-alphaIold)*np.dot(dataMatrix[i:i+1,:], dataMatrix[i:i+1,:].T) - labelMat[j, 0]*(alphas[j, 0]-alphaJold)*np.dot(dataMatrix[i:i+1,:], dataMatrix[j:j+1,:].T)\n",
    "                b2 = b - Ej- labelMat[i, 0]*(alphas[i, 0]-alphaIold)*np.dot(dataMatrix[i:i+1,:], dataMatrix[j:j+1,:].T) - labelMat[j, 0]*(alphas[j, 0]-alphaJold)*np.dot(dataMatrix[j:j+1,:], dataMatrix[j:j+1,:].T)\n",
    "                if (0 < alphas[i, 0]) and (C > alphas[i, 0]): b = b1\n",
    "                elif (0 < alphas[j, 0]) and (C > alphas[j, 0]): b = b2\n",
    "                else: b = (b1 + b2)/2.0\n",
    "                #更新alpha对修改次数计数器\n",
    "                alphaPairsChanged += 1\n",
    "                #打印更新次数\n",
    "                #print(\"循环次数: {} alpha:{}, alpha对修改了 {} 次\".format(iter,i,alphaPairsChanged))\n",
    "        if (alphaPairsChanged == 0): iter += 1\n",
    "        else: iter = 0\n",
    "        #print(\"迭代次数: {}\".format(iter))\n",
    "    return b,alphas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b = [[-3.84505534]]\n",
      "大于0的alpha\n",
      "[ 0.1277133   0.24197161  0.36968491]\n"
     ]
    }
   ],
   "source": [
    "dataArr, labelArr = loadDataSet(\"testSet.txt\")\n",
    "b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)\n",
    "print(\"b = {}\".format(b))\n",
    "print(\"大于0的alpha\\n{}\".format(alphas[alphas > 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.658191, 3.507396] -1.0\n",
      "[3.457096, -0.082216] -1.0\n",
      "[6.080573, 0.418886] 1.0\n"
     ]
    }
   ],
   "source": [
    "for i in range(alphas.shape[0]):\n",
    "    if alphas[i, 0] > 0.0:\n",
    "        print(dataArr[i], labelArr[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1.3.可视化数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualizeData(dataArr, labelArr, b, alphas):\n",
    "    \"\"\"\n",
    "    可视化数据\n",
    "    参数：\n",
    "        dataArr -- 参数矩阵\n",
    "        labelArr -- 标签矩阵\n",
    "        b -- 超平面的常数项\n",
    "        alphas -- 超平面的参数\n",
    "    返回：\n",
    "        无 -- 直接绘图\n",
    "    \"\"\"\n",
    "    dataMat = np.array(dataArr)\n",
    "    labelMat = np.array(labelArr).reshape(-1, 1)\n",
    "    posIndex = np.squeeze(labelMat > 0.0)\n",
    "    negIndex = np.squeeze(labelMat < 0.0)\n",
    "    svIndex = np.squeeze(alphas > 0.0)\n",
    "    posData = dataMat[posIndex, :]\n",
    "    negData = dataMat[negIndex, :]\n",
    "    svData = dataMat[svIndex, :]\n",
    "    plt.scatter(posData[:, 0], posData[:, 1], marker = 's', color = 'b')\n",
    "    plt.scatter(negData[:, 0], negData[:, 1], marker = '.', color = 'r')\n",
    "    plt.scatter(svData[:, 0], svData[:, 1], s = 200, marker = 'o', c='',edgecolors='g')\n",
    "    plt.legend((\"postive point\", \"negative point\", \"support vector\"))\n",
    "    w = np.dot((alphas * labelMat).T, dataMat)\n",
    "    x0 = np.array([2, 8])\n",
    "    x1 = -(w[0, 0] * x0 + np.squeeze(b)) / w[0, 1]\n",
    "    plt.plot(x0, x1, color = 'y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOXVwPHfk0lCWMIWdggEkC07EClRERC0sihqQaFo\nXaqIr7Ral1e0Lqhvq1WrtEWlrqCioGgVBa2KIqJYCGizsUOAsIYQskBCtvP+MUkMIcskM5M7y/n6\nmc/M3Ln3PmeGeM+9z3aNiKCUUsr/BFgdgFJKKWtoAlBKKT+lCUAppfyUJgCllPJTmgCUUspPaQJQ\nSik/pQlAKaX8lCYApZTyU5oAlFLKTwVaHUB9OnXqJBEREVaHoZRSXmPTpk3HRKSzI+t6dAKIiIgg\nKSnJ6jCUUsprGGP2OrquVgEppZSf0gSglFJ+ShOAUkr5KYfbAIwxrwGTgaMiEl2xrCOwDIgAMoCr\nRSSnlm0vBf4G2IBXROTJpgZcUlJCZmYmRUVFTd2FaiYhISH06tWLoKAgq0NRStWiMY3Ai4AFwBvV\nls0FVovIk8aYuRXv76u+kTHGBjwPXAxkAhuNMStEJL0pAWdmZhIaGkpERATGmKbsQjUDESE7O5vM\nzEz69u1rdThKqVo4XAUkImuB4zUWTwEWV7xeDFxRy6YjgJ0isltEioGlFds1SVFREWFhYXrw93DG\nGMLCwvRKTSkP5mwbQFcROVTx+jDQtZZ1egL7q73PrFhWK2PMLGNMkjEmKSsrq651mhiuak7676SU\nZ3NZI7DY7y3p9P0lReQlEUkQkYTOnR0ay6CU8lN79/6J/PzNVofhtZxNAEeMMd0BKp6P1rLOASC8\n2vteFcv82ocffkh6+s/NIA8//DBffvlls8exYsUKnnyy/jb5jIwM3n777WaKSCnHHDjwPHv2PMjR\no+9YHYrXcjYBrACur3h9PfBRLetsBAYYY/oaY4KB6RXb+bWaCeCxxx5j/PjxzR7H5Zdfzty5c+td\nRxOA8jTHjn3Mjh2/Jyzscvr1a3KnQr/ncAIwxrwDrAcGGWMyjTG/BZ4ELjbG7ADGV7zHGNPDGLMK\nQERKgTnAv4EtwLsikubar1G7tm3BmLMfbds2fZ8ZGRkMHjyYmTNnMmTIEKZOncqpU6cAWL16NUOH\nDiUmJoabbrqJ06dPAzB37lwiIyOJjY3lnnvu4fvvv2fFihXce++9xMfHs2vXLm644QaWL1/OZ599\nxrRp06rKW7NmDZMnTwbg888/JzExkWHDhjFt2jQKCgrOim/MmDHccccdxMfHEx0dzYYNGwA4fvw4\nV1xxBbGxsYwcOZLk5GQAFi1axJw5cwC44YYb+P3vf895551Hv379WL58eVX83377LfHx8Tz33HNN\n//GUcoG8vCTS06cTGjqMyMi3sXc0VE0iIh77GD58uNSUnp5+1rK6QN2PptqzZ48Asm7dOhERufHG\nG+Xpp5+WwsJC6dWrl2zbtk1ERK677jp57rnn5NixYzJw4EApLy8XEZGcnBwREbn++uvlvffeq9pv\n5fuSkhIJDw+XgoICERGZPXu2vPnmm5KVlSWjRo2qWv7kk0/Ko48+elZ8o0ePlptvvllERL755huJ\niooSEZE5c+bIvHnzRERk9erVEhcXJyIir7/+utx+++1VMUydOlXKysokLS1N+vfvLyIiX3/9tUya\nNKlJv1dj/r2UakhhYYasW9dV1q+PkNOnD1sdjkcCksTBY6yOBG6C8PBwzj//fACuvfZa1q1bx7Zt\n2+jbty8DBw4E4Prrr2ft2rW0a9eOkJAQfvvb3/LBBx/QqlWrevcdGBjIpZdeyscff0xpaSkrV65k\nypQp/PDDD6Snp3P++ecTHx/P4sWL2bu39jmfZsyYAcCFF15IXl4eJ06cYN26dVx33XUAXHTRRWRn\nZ5OXl3fWtldccQUBAQFERkZy5MiRJv9GSrlaSckJkpMnInKamJhVBAfX1ulQNYZHzwbqqWp2b6yv\nu2NgYCAbNmxg9erVLF++nAULFvDVV1/Vu//p06ezYMECOnbsSEJCAqGhoYgIF198Me+803CDV2Pi\nq6lFixZVr+0nE0pZr7y8mLS0qygs3EFs7Oe0bj3E6pB8gl4BNMG+fftYv349AG+//TYXXHABgwYN\nIiMjg507dwLw5ptvMnr0aAoKCsjNzWXixIk899xz/Pe//wUgNDSU/Pz8Wvc/evRoNm/ezMsvv8z0\n6dMBGDlyJN99913V/k+ePMn27dtr3X7ZsmUArFu3jnbt2tGuXTtGjRrFkiVLAHu7QqdOnWjrYGNI\nfbEq5W4iwrZtN3PixNcMGvQaHTqMsTokn6EJoAkGDRrE888/z5AhQ8jJyeG2224jJCSE119/nWnT\nphETE0NAQACzZ88mPz+fyZMnExsbywUXXMCzzz4L2M/yn376aYYOHcquXbvO2L/NZmPy5Ml8+umn\nVQ3AnTt3ZtGiRcyYMYPY2FgSExPZunVrrfGFhIQwdOhQZs+ezauvvgrAvHnz2LRpE7GxscydO5fF\nixfXum1tYmNjsdlsxMXFaSOwanYZGY9y5MibREQ8Trdu11odjk8xnnyZn5CQIDVvCLNlyxaGDHHs\n8q9tW6jtxDU0FGqp/nZIRkYGkydPJjU1tWk7cLMxY8bwzDPPkJCQYHUoQOP+vZSq6dChRWzbdiPd\nut3EoEGv6OhyBxhjNomIQwcAn24DaOpBXillvZyc1WzffgsdOoxn4MCFevB3A59OAO4QERHhsWf/\nYK/fV8rbFRSkkpp6Fa1aDSEqajkBATqluDtoG4BSyqOcPn2QlJSJ2GytiYlZSWBgO6tD8ll6BaCU\n8hilpQWkpFxGaWkO8fFrCQkJb3gj1WSaAJRSHqG8vJT09OkUFPyXmJiPCQ0danVIPk8TgFLKciLC\nzp2/5/jxlQwcuJCwsAlWh+QXtA3AA504cYIXXnih6v3BgweZOnWqJbFMnDiREydO1LvOokWLOHjw\nYDNFpHzR/v1/5eDBFwkPv48ePW61Ohy/oQnAA9VMAD169KiambO5rVq1ivbt29e7jiYA5YyjR99j\n9+576dz5Gvr1+7PV4fgVTQCNlJGRwZAhQ7jllluIiorikksuobCwEIBdu3Zx6aWXMnz4cEaNGlU1\nUnfXrl2MHDmSmJgYHnzwQdq0aQNAQUEB48aNY9iwYcTExPDRR/bbKcydO5ddu3YRHx/PvffeS0ZG\nBtHR0YB9Soi0tJ9n0x4zZgxJSUmcPHmSm266iREjRjB06NCqfVW3Zs0aLrzwQiZNmsSgQYOYPXs2\n5eXlALzzzjvExMQQHR3NfffdV7VNREQEx44dq/N7L1++nKSkJGbOnEl8fHzVb6GUI3Jzv2fLluto\n2/Z8Bg9ehDF6SGpWjk4basXD2emgq3z/vcif/2x/dtKePXvEZrPJjz/+KCIi06ZNkzfffFNERC66\n6CLZvn27iIj88MMPMnbsWBERmTRpkrz99tsiIvLiiy9K69atRUSkpKREcnNzRUQkKytL+vfvL+Xl\n5bJnz56qaZwry6x8/+yzz8rDDz8sIiIHDx6UgQMHiojI/fffXxVHTk6ODBgwoGrq6Epff/21tGjR\nQnbt2iWlpaUyfvx4ee+99+TAgQMSHh4uR48elZKSEhk7dqz861//EhGRPn36SFZWVr3fe/To0bJx\n48Zafy+dDlrV5eTJ7fLtt2Hyww8DpLj4mNXh+AwaMR207zcCr18P48ZBcTEEB8Pq1ZCY6NQu+/bt\nS3x8PADDhw8nIyODgoICvv/++zNu5lJ5Q5j169fz4YcfAvDrX/+ae+65B7An3wceeIC1a9cSEBDA\ngQMHGpyC+eqrr+aSSy7h0Ucf5d13361qG/j8889ZsWIFzzzzDABFRUXs27fvrGkYRowYQb9+/QD7\ntNHr1q0jKCiIMWPGUHkP5pkzZ7J27VquuOKKBr+3Uk1RXHyMlJSJGGOIiVlFUFCY1SH5Jd9PAGvW\n2A/+ZWX25zVrnE4A1adMttlsFBYWUl5eTvv27fnpp58c3s+SJUvIyspi06ZNBAUFERERQVFRUb3b\n9OzZk7CwMJKTk1m2bBkLFy4E7Mnk/fffZ9CgQfVu76qpoiu/t1KNVVZWSGrqFE6fziQu7itatTrH\n6pD8lu9XuI0ZYz/zt9nsz2PGuKWYtm3b0rdvX9577z3AfkCunPp55MiRvP/++wAsXbq0apvc3Fy6\ndOlCUFAQX3/9ddUNXhqafvmaa67hqaeeIjc3l9jYWAB++ctf8o9//KNqDv8ff/yx1m03bNjAnj17\nKC8vZ9myZVxwwQWMGDGCb775hmPHjlFWVsY777zD6NGjHf7uOl20cpRIOVu3Xk9e3nqGDHmLdu2c\nOxlTznE6ARhjBhljfqr2yDPG3FljnTHGmNxq6zzsbLkOS0y0V/s8/rhLqn/qs2TJEl599VXi4uKI\nioqqaoidP38+zz77LLGxsezcuZN27exD22fOnElSUhIxMTG88cYbDB48GICwsDDOP/98oqOjuffe\ne88qZ+rUqSxdupSrr766atlDDz1ESUkJsbGxREVF8dBDD9Ua47nnnsucOXMYMmQIffv25corr6R7\n9+48+eSTjB07lri4OIYPH86UKVMc/t433HADs2fP1kZg1aDdu+eSlfUe/fs/TefOv7I6HOVoY4Ej\nD8AGHAb61Fg+BviksftzWSOwxU6ePFl1T+B33nlHLr/8ckvicObevk3ljf9eyj0yM1+Qr79Gtm27\nver/B+V6WNgIPA7YJSK136zWT23atIk5c+YgIrRv357XXnvN6pCUalbZ2SvZsWMOYWGTOeec+Tq1\ns4dwdQKYDtR109rzjDHJwAHgHhFJq2M9nzNq1Kiq9gArjRkzhjFuagNRqi75+ZtIS7uGNm2GEhm5\nlIAA3+974i1c1ghsjAkGLgfeq+XjzUBvEYkF/gF8WM9+ZhljkowxSVlZWa4KTyllgaKifaSkTCYo\nKIyYmI+x2VpbHZKqxpWpeAKwWUTO6sguInnVXq8yxrxgjOkkIsdqWfcl4CWw3xKyqcFk5mXy8qaX\nWb1nNbmnc2kd1Jpf9PwFt517G4M7DW7qbpVSDiopOUFy8kTKygoZNuxLWrTobnVIqgZXdgOdQR3V\nP8aYbqai0s8YM6Ki3GwXll0l73QeMz+YSeyLsRwvPM7/XfR/vHXlW/z1kr8S2iKUMYvGcMmbl5CZ\nl+mO4pVSQHl5MWlpv6KwcDvR0R/QunWU1SGpWrjkCsAY0xq4GLi12rLZACKyEJgK3GaMKQUKgekV\nrdUudaLoBGMXj+UXPX/Bvj/so01wmzM+P7/3+Tw8+mGe+f4Zznv1PL654Rv6dujr6jCU8msiwrZt\nszhx4isGD15Mhw4XWR2SqoNLrgBE5KSIhIlIbrVlCysO/ojIAhGJEpE4ERkpIt+7otyaZn4wk1G9\nR/HipBfPOvhXCrYF88CoB7j3vHuZ9PYkSspK3BFKs5g/fz6nTp1y2/5/+uknVq1a5bb9K9+0d+/j\nHDmymIiIeXTr9hurw1H18JmRwD8d/omUIyk8+8tnHepi9rtf/I7OrTvz0bazZ830BmVlZR6ZAEpL\nS90UjfIGhw+/QUbGI3Ttej19+jTfeE/VND6TAF7c+CK3Dr+VwEZ0Mbv93Nt5fuPzjSrn5MmTTJo0\nibi4OKKjo1m2bBnw87TJAElJSVXdLefNm8d1111HYmIiAwYM4OWXXwaaNjVzmzZtuPvuu4mLi+NP\nf/oTBw8eZOzYsYwdO/aMGD/77LMzJqVbs2YNkydPBuyTxiUmJjJs2DCmTZtGQUEBABs3buS8884j\nLi6OESNGkJuby8MPP8yyZcuIj49n2bJlHD9+nCuuuILY2FhGjhxJcnLyGd/x/PPP57rrrmvU76l8\nR07OV2zb9lvat7+IQYNe0r7+3sDREWNWPBozErjr010lIyejrsFxtTpdelpa/6m1nCg84fA2y5cv\nl5tvvrnq/YkT9m0rp00WEdm4caOMHj1aREQeeeQRiY2NlVOnTklWVpb06tVLDhw40KSpmQFZtmxZ\nVdnVy6yupKREwsPDq6aDnj17trz55puSlZUlo0aNqlr+5JNPyqOPPiqnT5+Wvn37yoYNG0REJDc3\nV0pKSuT111+X22+/vWq/c+bMkXnz5omIyOrVqyUuLq7qOw4bNkxOnTp1Viw6Etg/FBSkytq17eQ/\n/4mS4uIcq8PxazRiJLDPXAFkF2bTPbRx3cyCbcGEtQojpyjH4W1iYmL44osvuO+++/j222+r5vWp\nz5QpU2jZsiWdOnVi7NixbNiwAfh5amabzVY1NfPGjRurpmYODAysmpoZ7DNw/upXDc+fEhgYyKWX\nXsrHH39MaWkpK1euZMqUKfzwww+kp6dz/vnnEx8fz+LFi9m7dy/btm2je/funHvuuYB9YrvAwLOv\npNatW1d1hn/RRReRnZ1NXp69h+/ll19Oy5YtHfsRlU85ffowyckTsdlaEhu7kqCg+u8gpzyHzwzJ\nCwkMoai0iGBbcKO2KyotIiQwxOH1Bw4cyObNm1m1ahUPPvgg48aN4+GHHyYwMLCqCqfmlM51TcHc\n2KmZQ0JCsNlsDsU5ffp0FixYQMeOHUlISCA0NBQR4eKLL+add87srZuSkuLQPuvTurUO8PFHpaUF\npKRMpqQkm6FD1xIS0sfqkFQj+MwVQFTnKL7d+22jttmevR2ATq06ObzNwYMHadWqFddeey333nsv\nmzdvBuxtAJs2bQKomvq50kcffURRURHZ2dmsWbOm6kzb2amZ65uGefTo0WzevJmXX36Z6dOnA/Zp\nqb/77jt27twJ2Nsztm/fzqBBgzh06BAbN24EID8/n9LS0rP2P2rUKJYsWQLY2xU6depE27ZtHf7t\nlG8RKWPLlhkUFPxIZORSQkOHWR2SaiSfSQC3Dr+VF5NebNQ2L258kZvib2pUw3FKSgojRowgPj6e\nRx99lAcffBCARx55hDvuuIOEhISzztJjY2MZO3YsI0eO5KGHHqJHjx6A81Mzz5o1i0svvfSsRmCw\nVxdNnjyZTz/9tKoBuHPnzixatIgZM2YQGxtLYmIiW7duJTg4mGXLlvG73/2OuLg4Lr74YoqKihg7\ndizp6elVjcDz5s1j06ZNxMbGMnfuXBYvXuzw76Z8i4iwY8cdZGd/woAB/6BTp8lWh6SawtHGAise\njWkEPll8Ujo91Uk2ZG5wqKFk74m90vEvHWVPzh6H1m+qRx55RJ5++umzllsxNbMVtBHYN+3b91f5\n+mtk5857rA5F1YA/NgK3CmrFK5e9wpSlU0g+klzvugfyDnDpW5fywAUPENE+onkCVMpHZGW9z65d\n99C58zT69fuL1eEoJ/hMIzDAlMFTKCotYvwb47l52M3cOvxW+rT/uVEq62QWr/34Gn/7z9+4c+Sd\n3JV4l9tjmjdvXq3LdWpm5Y1yc9ezZcu1tG2byODBizHGZ84h/ZJXJgARqbPHzDXR1zCs+zBe2PgC\nw14axjkdz6Fjy47kn84n9WgqVw25io9nfMzwHsObOWr/I66f7klZ6NSpnaSmXk6LFr2Ijv4Im027\n/Xo748n/kyYkJEhSUtIZy/bs2UNoaChhYWENdps8WXySHw//SN7pPFoHtSamawwdW3Z0Z8iqgoiQ\nnZ1Nfn4+ffvqhHverqQkm82bEykpOc6wYetp1WqA1SGpOhhjNolIgiPret0VQK9evcjMzMTRm8WE\nVfxHKRzJOMIRzrpdgXKTkJAQevXqZXUYykllZUWkpEyhqGgf8fFf6cHfh3hdAggKCtIzSqWaiUg5\nW7feQF7ed0RGvku7dudZHZJyIW3BUUrVaffuB8jKWka/fk/Rpcu0hjdQXkUTgFKqVgcP/pP9+/9C\njx63ER5+j9XhKDfQBKCUOkt29qds3347HTtO4pxz/q5TO/soTQBKqTPk5/9IWto02rSJJTJyKQGN\nmCpFeRdNAEqpKkVF+0lJmUxQUEdiYj4hMLD2W6sq3+CSBGCMyTDGpBhjfjLGJNXyuTHG/N0Ys9MY\nk2yM0WkDlfIwpaW5pKRMpKysgJiYVbRo0cPqkJSbufLabqyIHKvjswnAgIrHL4AXK56VUh6gvLyE\ntLRpnDq1ldjYz2jTJtrqkFQzaK4qoCnAGxWT1f0AtDfGNO72XUoptxARtm+fTU7OFwwc+DIdOoyz\nOiTVTFyVAAT40hizyRgzq5bPewL7q73PrFimlLLY3r1/4vDh1+jT52G6d7/B6nBUM3JVFdAFInLA\nGNMF+MIYs1VE1jZlRxUJZBZA7969XRSeUqo2hw+/RUbGQ3Tteh0REfOsDkc1M5dcAYjIgYrno8C/\ngBE1VjkAhFd736tiWW37eklEEkQkoXPnzq4ITylVi5ycNWzbdhPt249l0KBXtK+/H3I6ARhjWhtj\nQitfA5cAqTVWWwH8pqI30EggV0QOOVu2UqppTp7cQlralbRsOYCoqA8ICAi2OiRlAVdUAXUF/lVx\n9hAIvC0inxljZgOIyEJgFTAR2AmcAm50QblKqSYoLj5CSspEjGlBTMxKgoLaWx2SsojTCUBEdgNx\ntSxfWO21ALc7W5ZSyjllZSdJSZlMcfFR4uO/oWXLCKtDUhbSMd5K+QmRMtLTf01+/maioz+kbVuH\n7hmifJgmAKX8xM6dd5GdvYIBAxbQqdNlVoejPIDOBaSUH9i/fz4HDvydXr3uomdPrY1VdnoFoFQT\nHC88ztLUpezO2U1ZeRnd2nTj6qir6dvB8+5Wl5X1L3btuotOnX5F//5PWx2O8iB6BaC82/r18MQT\n9udmsC93Hzd9dBP9/96f7/Z/R9fWXQlvF87+vP2c+/K5THp7Ev/J/E+zxOKIvLz/sGXLr2nb9hcM\nGfImxuj/8upnegWgvNf69TBuHBQXQ3AwrF4NiYluKy75SDITlkzgxvgb2TZnG11adznj86cufool\nyUu47J3LeHHSi/wq8ldui8URhYW7SUm5jODgHkRHr8Bma2lpPMrzaAJQ3mvNGvvBv6zM/rxmjdsS\nwMH8g0xcMpFnLn6GGTEzal2nVVArbhl+Cwk9EvjlW7+ka5uuXND7ArfE05CSkuMkJ09EpIzY2E8J\nDtZR9epsej2ovNeYMfYzf5vN/jxmjNuKenb9s0yNnFrnwb+6od2HMv/S+Tyw+gG3xVOf8vLTpKZe\nQVHRHqKjP6JVq4GWxKE8nyYA5b0SE+3VPo8/7tbqn8KSQhb/dzG//8XvHd5mWuQ0duXsIuVIilti\nqotIOVu33khu7rcMGfIG7dtbcwWivIMmAOXdEhPh/vvdWve/csdKhnUfRr8O/c7+sI5G6CBbEDfF\n38SbyW+6La7a7NnzEEePvkO/fk/Spcs1zVq28j7aBqBUAzLzMhkcNvjsDxpohB7UaRCf7vy02eI8\nePBl9u37M927zyI8/H+brVzlvfQKQCkHCHL2wtoaoatvI7Vs4ybHj/+b7dtvo2PHSxkw4Hmd2lk5\nRBOAUg0IbxvOlmNbzv6ggUborce20iu0l9vjKyj4L2lpU2nTJobIyHcJCNALe+UYTQBKNWDSwEn8\ndPgndh3fdeYH9TRCF5cV89pPr/GbuN+4NbaiokySkycRGNiemJhPCAwMdWt5/qBtWzDm7EfbtlZH\n5nqaAJRqQEhgCDfG38j8H+af/WEdjdDLUpcxKGwQUV2i3BZXaWkeKSmTKCvLIyZmJS1a6G22XSE/\nv3HLvZkmAKUccFfiXXy07SPe+O8bDa6bdDCJuz+/myfGPeG2eMrLS0hLm8apU+lERb1PmzaxbitL\n+S5NAEo5oFubbnw681P++NUfue+L+ziUf/YdTQuKC3hh4wtMWDKBVy5/hcRw93RNFRF27PgfcnI+\nZ+DAf9Kx48VuKcdX+VMVT0O0tUgpB0V1ieKH3/7A42sfJ/KFSMb1HUd8t3hsxsbunN28v+V9RkeM\n5rOZnzG8x3C3xbFv3xMcOvQKffo8SPfuN7mtHF/lT1U8DTHOdlUzxoQDb2C/N7AAL4nI32qsMwb4\nCNhTsegDEXmsoX0nJCRIUlKSU/Ep5Q65Rbm8m/YuGScyKC0vpVubbkyNnEp4u3C3lnvkyNts2TKT\nLl1mVszuqd09G6u+n0yk4c89nTFmk4g4dLs3V1wBlAJ3i8hmY0wosMkY84WIpNdY71sRmeyC8pSv\nWb/e3od+zBi3juh1pXYh7bhl+C3NWuaJE2vZuvVG2rUbzeDBr+rB301CQ2u/Ggj1wQ5Wrrgp/CHg\nUMXrfGPMFqAnUDMBKHW2Zp7S2VudPLmV1NQraNmyH9HR/yIgoIXVIfmsvDyrI2g+Lm0ENsZEAEOB\n2u6IcZ4xJtkY86kxxn1945R3aWA0rYLi4iOkpEzEmCBiYlYRFNTB6pCUj3BZAjDGtAHeB+4UkZo5\ndDPQW0RigX8AH9azn1nGmCRjTFJWVparwlOeqhmndPZGZWWnSEm5nOLiw8TEfELLlp53y0lvU1dV\nji9W8TTE6UZgAGNMEPAJ8G8RedaB9TOABBE5Vt962gjsJ7ywDaA5iJSRljaVY8c+Ijr6X3TqNMXq\nkJQXaNZGYGNviXoV2FLXwd8Y0w04IiJijBmB/coj29mylY9ITHTNgd/HEsmuXfdw7NiHnHPO3/Tg\n7wZt29bd2NuYdgBX7ccKrugFdD5wHZBijPmpYtkDQG8AEVkITAVuM8aUAoXAdGnOqRKV7/OxxuTM\nzL+TmTmfXr3upFcvx29EoxznqvEA3jyuwBW9gNYB9fZHE5EFwAJny1LqLJVn/fv2Ndv9gd3t2LGP\n2LnzTjp1upL+/Z+xOhzlw3QksPJe1c/6bTYIrPhz9uLG5Ly8DaSnzyA09FyGDHkLY2xWh6R8mCYA\n5b2qdyEFuOUW6N3ba9sACgv3kJJyGcHB3YiJ+RibrZXVISkfp5PBKdeq4x65blGzC+nQoe4v001K\nSnJISZmISAmxsZ8SHNzF6pB8Rl2Tvym9AlDOqNnrprkbYitvyLJmDYSFwZ13emUjcHn5aVJTr6Sw\ncDdxcV/QqtUgq0PySnX1xmmsxo4HcHbqCCt7EWkCUE1T28G+tlG97j4IV3YhfeIJr2wEFhG2bv0t\nubnfMGTiEmTlAAAczUlEQVTI27Rvf6HVIXktZw/+Te2X6OxB2speRFoFpGrXUFXOmjVw+rT9gHv6\n9M9XAlaN6vXSEcUZGQ9z9OgS+vb9M127zrA6HOUEd9xnwN33KdArAHU2R6pywsKgvNz+urzc/r56\nlYyrGmIdHdzljrLd7NCh19i79//o3v1meveea3U4XstVVT/OcteZvDu/myYAdTZHqnKysyEgwH7w\nDwiwvwfXjeqFxrcpuLJsNzt+/HO2bZtFhw6/ZMCAF3RqZyd4wsHfW2kVkDqbI9UpY8ZAixb2dVq0\ncE+Vi4/OFFpQkExa2lRat44iKupdAgKCrA5J+Sm9AlBnc6Q6pTmqXCoTUeUVgJfU69fn9OkDpKRM\nwmZrS0zMSgID/fBGtM2sem+a+nrc1MedPXXq6kXUHFwyG6i76GygypcmeCstzefHH0dRVLSboUO/\npU2bOLeWl5mXyQdbPuDoyaMEBQQxIGwAVw6+kpZBLd1abnNwtN7fVYc3R24T6eytJF11K8rGzAaq\nVUDK9Vw5GCwxEe6/3+sP/uXlpaSnX83Jk6lERb3n1oP/xgMbuWrZVcQtjCP5SDLBtmDKpIy3kt+i\n9/ze3PP5PWSf8u7JeB09Y3ZnD5qanL3PgBX3KdAqIOW86mfp4NxgMB86468kIuzYcTvHj3/GwIEv\n07HjL91W1vL05fzPyv/h0TGP8saVb9AmuM0Zn+/O2c1T3z1F4quJfH7d50S0j3BbLJ4gP//nM2t3\nD6xydt9WTB2tCUA5p+aEbPHx9nEB5eWNH5DlY1M6V9q//ykOHXqJ3r0foEePm91WzpqMNdy+6na+\nuO4L4rrVfoXRr0M/Fk5eyHPrn2PCkgms/+162oe0d1tMnkR7C51Nq4D8javn6nnjDSgq+rmnzoYN\nP3cNbWzDrQ/2+jlyZCm7d8+lS5cZ9O37uNvKERHu/vxuFk5aWOfBv7o/JP6BYd2HsTBpodtiUp5P\nE4A/qTzDfugh+3NlEmhqUli/Hl5//ewWqoAAGD++8WfwjRjN68mdFyqdOLGOrVuvp127UQwe/DrG\nuO9/t40HN5JTmMOUwT/fOayhkal3jbyLhUkLKSsvc3gbK1WPrbn56n2EtQrIn9R1ht3Uapc1a6C0\n1P7aGPuBH+z7mTev8dU39XQtFRE2HtzICxtf4OPtH5NTmENIYAjRXaK5dfitzIiZQasgz5k++dSp\nbaSmTiEkpC/R0R8SENDCreUt+mkRs4bPIqBakmloZOrwHsPp1KoTazLWMK7fOIe2sZKVMeTl1d7z\nKD/fvtzTb/1YF70C8Ce1nWE7U+1SfX8hIfDCC/D44y6vu88+lc34N8dzzfJriOwcScptKZQ8VELW\nvVk8OuZRPtz2Ib2f682KbStcVqYziouzSE6eiDE2YmNXERTU0e1l7s/bT1TnqEZvF9k5ksy8TDdE\n9LOGrio8+aqjOk9Ojk2lVwDeqCk9ZSq3mT/fPm1D9W0bO9iqevmuHAxWSyNwTvxgLlx0IZMGTOLz\naz/HFvDzHbJaB7dmwoAJTBgwgQ0HNnDF0is4XXqaaVHTnIvDCWVlhaSmXk5x8SHi47+mZct+zVJu\nuZQ3aTtjTJO3dVRDB05fPLB6C5ckAGPMpcDfABvwiog8WeNzU/H5ROAUcIOIbHZF2X6nKT1l6tum\nsSN6a9vX/fe74pvVejXy231Pc0m/S3jq4qfq3XREzxGsmrmKi9+8mGHdh9G/Y3/XxNQIIuVs2XIt\neXn/ISrqfdq2/UWzld0ztCc7ju9o9Hbbs7fz6+hfuyEiz1Sz/aA55tz3ZE5XARn7TUufByYAkcAM\nY0xkjdUmAAMqHrOAF50t1281pcqmoW1qDraqr1HYnT11alRR7R4xgLV71/KncX9yaPP4bvHcEHeD\nZT1bdu26l2PHPqB//2fp3PnKZi372threXnzy41qHE85ksK+3H1c1PciN0bm2fz9KsMVVwAjgJ0i\nshvAGLMUmAKkV1tnCvCG2P86fzDGtDfGdBeRQy4o3780ZX4cR7dZv97erfO11+wH+MBAuPFG+60W\nK6uN3DE/Tx1VSv8s+JDr465vVOPu7ITZjHx1JI+NfaxZpzzIzFxAZuaz9Oz5e8LD72y2ciuN6j0K\nm7Hx+a7P+eU59oFmDd2p6u//+Tu3DLuFIFvQGZ85c3erxtJJUC0mIk49gKnYq30q318HLKixzifA\nBdXerwYSGtr38OHDxe98/73In/9sf3ZmncZu8/33Ii1bihgjYu/YaX9Uvg8IsH/+/fdNK7++uFq2\nFLHZft5/hYSXEmTd3nWN/j5DFw6VH/b/4HxsDsrKWiFffx0gyclTpLy8tNnKrenjbR9L92e6y47s\nHQ2u+/qPr0uf5/rI0YKjTpUZGnrmn0vlIzT053Vq+7wxj8p9Obufuh6u/K6eAEgSB4/fHtcIbIyZ\nhb2aiN69e1scTTNztH6/tnnvG2oYrrlNzfUrq3ZqViFUvq8+stcVc/NUlr9vX533Hsg/nU+Hlh1q\n37ae36lDyw7knW6eit28vCTS06cTGjqcyMgl2GtErTF54GQeG/sYo14fxfxfzueqIVedcXYPcPTk\nUeb/MJ+3kt/i39f+m86tOztVZnM04DZ20rfGXlUY41hbgC+2FbgiARwAwqu971WxrLHrACAiLwEv\ngX02UBfE5z2aek/dl16C22+3H6RbtGi4Ybi2A2hYmL0fvwgEBcGECbBqlb2ff1NH9jpSfmCgvc4f\nztp/6+DWFBQXnL19A79TQXHBWXPguENhYQYpKZMJDu5CTMzH2Gyt3V5mQ24edjMR7SN4fO3j/OHf\nf+Da2Gvp3a43JWUlJB1KYuX2lVw55Eq+/+339Grby+pwG8Wd1VP+2hbgigSwERhgjOmL/aA+HajZ\nrWAFMKeifeAXQK5o/f/ZmlK/vn49zJnz84Csyvvz1pcAah5A33gDFi+278MYuOMO+Mtffj5LDws7\nu+uoM6qXD3DLLdC791n7H9lzJKt2rGJEzxFnbl/P73Qo/xA7sncQ1aXxfeIbo6TkBCkpExE5TUzM\n1wQHd3VreY0xvt94xvcbT+rRVJanLyc9K51gWzAje47kHxP+QceW7h+X4A6+eAZuNacTgIiUGmPm\nAP/G3g30NRFJM8bMrvh8IbAKexfQndi7gd7obLk+qSk3WVmz5ucDKdjP1OtLHOvX26tdqp91gz1x\nVFZrPvccXHGF+26xWPMA/pvf1FrObefexiVvXsIfR/3xzKqMen6nlze/zNVRV9O2hftGEZWXF5OW\ndhWFhTuJi/uC1q2HuK0sZ0R3iSa6S7TVYbidlfcEdueNYpqDS9oARGQV9oN89WULq70W4HZXlOV2\nVk9HXFlmZffKumKofnbeooX9AG6zwYIF9m1q+x41q15uucV+8AV45ZWfb/JeVlb/VYSzv5GDiS56\nVz6DTobw7LI7uO/XL5y9jxrbZZzI4PmNz/PFdV80PiYHiQjbtt3MiRNfM2TIW7RvP9ptZSnHWFl9\n4/WD2BxtLbbi0ey9gOrpkeJRMdRc55//PLNHTOXnAQEiQUH2z0Xs69hs9vN8m83+vtI//2lft3pv\nn6bG5woV5ezrECC97jKyYOnd9a6+M3unDPj7AJm/fr574qmwe/fD8vXXyJ49j7u1HG/hSM+Yutap\nbZum9LJpTG8iV/QGcrRsq9CIXkA6F1B1njAdsSMx1FwnO/vMnjlr1vw8J39Jib2BeP36+mfbnDUL\nvvkG/u//zmxErjkorLl+o4pywnPK+WaxYcHWN7jw9QtZmrqU4rLiqtVSj6Zy+8rbOfflc7lz5J3c\nMfIO98QDHDq0iL17H6Nbt5vo0+ePbivHm+Tl1X74q179UX0dZ/fVWDX346uzejaVx3UDtZQn3ITc\nkRgaWmfMGPtBvrJKp7z85+6b9VW91NZVtGZvoeb6jaqV0+9UMMnj3+ejDkd4YeML3PTRTXRo2YGi\n0iJa2Fpwy7BbSL4t2a29Wo4f/5Lt22+hQ4eLGThwIUZHMHmFml08vaFevjlpAqiuKY2wNTVH/XhD\n6yQm2tsCqncNrTxQO9rGULlOzbP9hpKIq9T4jkGJiUwFpkZOpaC4gJzCHFoGtaRDSIczJohzh4KC\nVNLSfkWrVkOIinqPgICghjfyUw01ijb3SGPwovp4Cxhp6LrMQgkJCZKUlGR1GI7ztFsaNtQQHBxc\n++yg1bf3pO9jgdOnD7J580hEyhg27AdCQsIb3siP1Xdh5K5DjSO9gJq7bCt7ARljNolIgiPr6hWA\nKzV1IJe71NaNs3qMp0/bxxCUl9d+gHfFFZEXKy0tICVlMqWlOcTHf6sHfw9VeaC1olbO26uUNAG4\nkie0IdSncgxAYMU/uzH2RFDfDdzdNRbAw5WXl5Kefg0FBcnExHxMaGi81SEp6j/jVo2nCcCVPPmM\nuXp1js1mHwMwdCjceafnJiyLiAg7d/6O48dXMXDgPwkLm2B1SKqC1/e79zCaAFzNU8+Ya06/0Lu3\nvetnTIxnJiwL7d//DAcPLiQ8/D569JhldTjKSXp1UDdNAP6iruopT01YFjl69F127/5fOne+hn79\n/mx1OF6nKb183NWQ6sH9WzyGJgB/4cnVUx4iN/c7tmz5De3aXcDgwYswRsdJNlZTDtharWMdTQC+\nwNGxB550tm/1nEs1nDq1g5SUKYSE9CY6+kNsthCrQ/IJrjy7b+u++f38lm8mAA87uLiVN/bV97CY\ni4uzSE6egDGG2NhPCQoKsywWX+PKs3u9InA937vGrTy4PPSQ/bm2G5v7ElfOzVPfzeBdyRPmXKpQ\nVlZIauoUiosPEB29gpYt+1sWi3Idbfh1jO9dAXjaYCx3c9XYg+Y8K/eQ8RIi5Wzd+hvy8n4gKuo9\n2rXz4b8TH6cNvk3jewnAQw4uzcZVjbvNmTg9pEF69+77yMpaTv/+f6Vz519ZEoOyZn4gZed7CcBD\nDi7NyhWNu82dOC1ukD5w4AX273+GHj1up1evP1gWh/L+6RS8me8lALD84OKV/ChxHjv2CTt2/I6w\nsMsYMOBvOrWzizkyOVtTzu71SsH1fDMBqKbxg8SZn7+J9PRraNNmKJGR72CMe6eSbognzibprPoO\n/s7U1Xvr7+HJnEoAxpingcuAYmAXcKOInKhlvQwgHygDSh2dqlQpVyoq2ktKymSCgjoTE/MJNltr\nq0Pyu0FQ1S+2vDnJ+Qpnu4F+AUSLSCywHbi/nnXHiki8HvyVFUpKTpCcPImyskJiY1fSokU3q0Py\ne76a5LyJUwlARD4XkdKKtz8A7rsnn1JNVF5eTFrarygs3E509Ae0bh1ldUhKeQRXDgS7Cfi0js8E\n+NIYs8kYU+/0isaYWcaYJGNMUlZWlgvDU/5IRNi2bRYnTnzFoEGv0KHDRVaHpJTHaLANwBjzJVDb\n9fIfReSjinX+CJQCS+rYzQUicsAY0wX4whizVUTW1raiiLwEvAT2W0I68B18iz9NY9EM9u59jCNH\nFhMR8Sjduv3G6nD8Ql29dZTnaTABiMj4+j43xtwATAbGSR03GBaRAxXPR40x/wJGALUmAL/mYXPk\neLvDhxeTkTGPbt1uoE+fh6wOp1a+2LWxZsOu9rL1XE5VARljLgX+F7hcRE7VsU5rY0xo5WvgEiDV\nmXJ9lgfNkePtcnK+Ytu2m2nffhwDB/7TY/v65+XZu0bWfPhS75i6kpk3Jzlf4WwbwAIgFHu1zk/G\nmIUAxpgexphVFet0BdYZY/4LbABWishnTpbrmypH49ps/jGNhZucPJlGaupVtGw5iOjo9wkICLY6\nJL/mD0nOWzk1DkBEzqlj+UFgYsXr3UCcM+X4DT8ajesup08fIjl5IjZbK2JjVxEY2M7qkJTyWDoS\n2NP4wWhcdyktLSAlZTIlJdkMHbqWkJDeVoeklEfTBKB8Qnl5KVu2zKCg4CdiYj4mNHSY1SEp5fE0\nASivJyLs3HkH2dmfMGDAi4SFTbQ6JKW8gu/dEUz5nczMZzl48AXCw++lZ8/ZVoejlNfQBKC82tGj\ny9m16x46d55Gv35PWh2OUl5FE4DyWrm569m69Tratj2PwYPfwBj9cwb7FNPGnP1o29bqyJSn0f9j\nlFc6dWonqamX06JFL6KjP8JmC7E6JI/hb1NMq6bTBKC8TnHxMVJSJiIixMR8SnBwJ6tDUsoraS8g\n5VXKyopITb2CoqJ9xMd/RatWtY5FVEo5QBOA8hoi5Wzdej15ed8RGfku7dqdZ3VISnk1rQJSXmP3\n7gfIynqXfv2epkuXaVaHo5TX0wSgvMLBg/9k//6/0KPHbYSH3211OB5NZ99UjtIqIOXxsrNXsX37\n/9Cx4yTOOefvHju1s6fQWTaVo/QKQHm0/PwfSUu7mjZt4omMXEpAgJ6zKOUqmgCUxyoq2k9KyiSC\ngjoSE/MJgYFtrA5JKZ+ip1PKI5WW5pKSMpGyspMMHfodLVp0tzokpXyOJgDlccrLS0hLm8qpU1uJ\njf2MNm2irQ5JKZ+kCUB5FBFh+/Zbycn5ksGDF9GhwzirQ1LKZzl7U/h5xpgDFfcD/skYU+tE7MaY\nS40x24wxO40xc50pU/m2vXv/xOHDr9OnzyN063a91eEo5dNccQXwnIg8U9eHxhgb8DxwMZAJbDTG\nrBCRdBeUrXzI4cNvkZHxEF27/oaIiEesDkcpn9ccvYBGADtFZLeIFANLgSnNUK7yIjk5a9i27Sba\ntx/LoEEva19/pZqBKxLA74wxycaY14wxHWr5vCewv9r7zIplSgFw8uQW0tKupGXLAURFfUBAQLDV\nISnlFxpMAMaYL40xqbU8pgAvAv2AeOAQ8FdnAzLGzDLGJBljkrKyspzdnfJwp08fJjl5AgEBIcTG\nriIoqL3VISnlNxpsAxCR8Y7syBjzMvBJLR8dAMKrve9Vsayu8l4CXgJISEgQR8pW3qms7CSpqZdR\nUpJFfPw3hIT0sTokpfyKs72Aqo/OuRJIrWW1jcAAY0xfY0wwMB1Y4Uy5yvuJlJGe/mvy8zcTGbmU\ntm0TrA5JKb/jbC+gp4wx8YAAGcCtAMaYHsArIjJRREqNMXOAfwM24DURSXOyXOXFRISdO/9AdvYK\nBgxYQKdOl1kdklJ+yakEICLX1bH8IDCx2vtVwCpnylK+IzPzbxw48A969bqbnj1vtzocpfyWTgan\nmlVW1gfs2nUXnTr9iv79n7I6HKX8miaAmtavhyeesD8rl8rL+w9btsykbdtfMGTImxijf35KWUnn\nAqpu/XoYNw6KiyE4GFavhsREq6PyCYWFu0hJuYzg4J5ER6/AZmtpdUhK+T09BatuzRr7wb+szP68\nZo3VEfmEkpLjJCdPRKSM2NhVBAd3tjokpRSaAM40Zoz9zN9msz+PGWN1RF6vrKyI1NQrKCraS3T0\nR7RqNdDqkJRSFbQKqLrERHu1z5o19oO/Vv84RaScbdtuJDf3WyIjl9K+/QVWh6SUqkYTQE2JiXrg\nd5E9ex7k6NGl9Ov3JF26XGN1OEqpGrQKSLnFwYMvsW/fE3Tvfivh4f9rdThKqVpoAlAul539Gdu3\n/w8dO05gwIAFOrWzUh5KE4Byqfz8n0hPn0abNjFERi4jIEBrGZXyVJoAlMsUFWWSkjKJwMD2xMSs\nJDAw1OqQlFL10NMz5RKlpXmkpEyirCyfoUO/o0WLHlaHpJRqgCYA5bTy8hLS0qZx6lQ6MTGraNMm\nxuqQlFIO0ASgnCIibN9+Gzk5nzNo0Kt07Hix1SEppRykbQDKKfv2PcHhw6/Sp8+DdO9+k9XhKKUa\nQROAarIjR95mz54/0rXrtUREPGZ1OEqpRtIEoJrkxIlv2Lr1Rtq3H8OgQa9oX3+lvJAmANVoJ09u\nJTX1Slq27EdU1AcEBLSwOiSlVBNoAlCNUlx8hJSUCRgTREzMKoKCOlgdklKqiZzqBWSMWQYMqnjb\nHjghIvG1rJcB5ANlQKmIJDhTrrJGWdkpUlIup7j4CPHx39CyZV+rQ1JKOcHZm8JXTfFojPkrkFvP\n6mNF5Jgz5SnriJSxZctM8vM3Eh39L9q2PdfqkJRSTnLJOABjbwG8GrjIFftTnmfnzrs5duxDzjnn\n73TqNMXqcJRSLuCqNoBRwBER2VHH5wJ8aYzZZIyZVd+OjDGzjDFJxpikrKwsF4WnnJGZ+TcOHPgb\nvXr9gV69fmd1OEopF2nwCsAY8yXQrZaP/igiH1W8ngG8U89uLhCRA8aYLsAXxpitIrK2thVF5CXg\nJYCEhARpKD7lXllZH7Jz5x/o1OlK+vd/2upwlFIu1GACEJHx9X1ujAkErgKG17OPAxXPR40x/wJG\nALUmAOU58vI2sGXLrwkNHcGQIW9hjM3qkJRSLuSKKqDxwFYRyaztQ2NMa2NMaOVr4BIg1QXlKjcq\nLNxDSsplBAd3JyZmBTZbK6tDUkq5mCsSwHRqVP8YY3oYY1ZVvO0KrDPG/BfYAKwUkc9cUK5yk5KS\n4yQnT0CkhNjYVQQHd7E6JKWUGzjdC0hEbqhl2UFgYsXr3UCcs+Wo5lFefprU1KsoKtpDXNyXtGo1\nqOGNlFJeSaeDVlVEhK1bbyI39xuGDHmb9u1HWR2SUsqNdCoIdYZWrYbQt++f6dp1htWhKKXcTK8A\nVBVjDBERD1odhlKqmegVgFJK+SlNAEop5ac0ASillJ/SBKCUUn5KE4BSSvkpTQBKKeWnNAEopZSf\n0gSglFJ+yoh47pT7xpgsYK/VcdTQCfD1W1vqd/Qd/vA99TueqY+IdHZkRY9OAJ7IGJPk6ze11+/o\nO/zhe+p3bDqtAlJKKT+lCUAppfyUJoDGe8nqAJqBfkff4Q/fU79jE2kbgFJK+Sm9AlBKKT+lCcBB\nxphLjTHbjDE7jTFzrY7HHYwx4caYr40x6caYNGPMHVbH5C7GGJsx5kdjzCdWx+IOxpj2xpjlxpit\nxpgtxphEq2NyNWPMHyr+TlONMe8YY0KsjskVjDGvGWOOGmNSqy3raIz5whizo+K5gyvK0gTgAGOM\nDXgemABEAjOMMZHWRuUWpcDdIhIJjARu99HvCXAHsMXqINzob8BnIjIY+z25feq7GmN6Ar8HEkQk\nGrAB062NymUWAZfWWDYXWC0iA4DVFe+dpgnAMSOAnSKyW0SKgaXAFItjcjkROSQimyte52M/aPS0\nNirXM8b0AiYBr1gdizsYY9oBFwKvAohIsYicsDYqtwgEWhpjAoFWwEGL43EJEVkLHK+xeAqwuOL1\nYuAKV5SlCcAxPYH91d5n4oMHxuqMMRHAUOA/1kbiFvOB/wXKrQ7ETfoCWcDrFdVcrxhjWlsdlCuJ\nyAHgGWAfcAjIFZHPrY3KrbqKyKGK14eBrq7YqSYAdRZjTBvgfeBOEcmzOh5XMsZMBo6KyCarY3Gj\nQGAY8KKIDAVO4qIqA09RUQc+BXuy6wG0NsZca21UzUPsXTdd0n1TE4BjDgDh1d73qljmc4wxQdgP\n/ktE5AOr43GD84HLjTEZ2KvyLjLGvGVtSC6XCWSKSOXV23LsCcGXjAf2iEiWiJQAHwDnWRyTOx0x\nxnQHqHg+6oqdagJwzEZggDGmrzEmGHtj0wqLY3I5Y4zBXm+8RUSetToedxCR+0Wkl4hEYP93/EpE\nfOrMUUQOA/uNMYMqFo0D0i0MyR32ASONMa0q/m7H4WMN3TWsAK6veH098JErdhroip34OhEpNcbM\nAf6NvbfBayKSZnFY7nA+cB2QYoz5qWLZAyKyysKYVNP8DlhSccKyG7jR4nhcSkT+Y4xZDmzG3nvt\nR3xkRLAx5h1gDNDJGJMJPAI8CbxrjPkt9hmSr3ZJWToSWCml/JNWASmllJ/SBKCUUn5KE4BSSvkp\nTQBKKeWnNAEopZSf0gSglFJ+ShOAUkr5KU0ASinlp/4ffDfU4MFGl6wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1aecb5fc7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualizeData(dataArr, labelArr, b, alphas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2.完整版Platt SMO的支持函数  \n",
    "\n",
    "### 6.2.1.Platt SMO支持函数  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class optStructK:\n",
    "    def __init__(self,dataMatIn, classLabels, C, toler):\n",
    "        self.X = dataMatIn\n",
    "        self.labelMat = classLabels\n",
    "        self.C = C\n",
    "        self.tol = toler\n",
    "        self.m = np.shape(dataMatIn)[0]\n",
    "        self.alphas = np.zeros((self.m,1))\n",
    "        self.b = 0\n",
    "        #误差缓存\n",
    "        self.eCache = np.zeros((self.m,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calcEkK(oS, k):\n",
    "    fXk = np.float(np.dot((oS.alphas * oS.labelMat).T, np.dot(oS.X, oS.X[k,:].T))) + oS.b\n",
    "    Ek = fXk - np.float(oS.labelMat[k])\n",
    "    return Ek"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def selectJK(i, oS, Ei):\n",
    "    \"\"\"\n",
    "    内循环选择alpha的启发式算法\n",
    "    参数：\n",
    "        i -- 外循环alpha的下标\n",
    "        oS -- 类\n",
    "        Ei -- 误差\n",
    "    返回：\n",
    "        j -- 选择的alpha的下标\n",
    "        Ej -- 误差\n",
    "    \"\"\"\n",
    "    #初始化\n",
    "    maxK = -1; maxDeltaE = 0; Ej = 0\n",
    "    oS.eCache[i] = [1,Ei]\n",
    "    #选择合理的集合\n",
    "    validEcacheList = np.nonzero(oS.eCache[:,0])[0]\n",
    "    if (len(validEcacheList)) > 1:\n",
    "        #选择最大步长的alpha\n",
    "        for k in validEcacheList:\n",
    "            #不重复计算\n",
    "            if k == i: continue\n",
    "            #计算误差\n",
    "            Ek = calcEkK(oS, k)\n",
    "            #计算步长\n",
    "            deltaE = abs(Ei - Ek)\n",
    "            #记录最佳选择\n",
    "            if (deltaE > maxDeltaE):\n",
    "                maxK = k; maxDeltaE = deltaE; Ej = Ek\n",
    "        return maxK, Ej\n",
    "    #没有合理值\n",
    "    else:\n",
    "        #随机选择\n",
    "        j = selectJrand(i, oS.m)\n",
    "        Ej = calcEkK(oS, j)\n",
    "    return j, Ej"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def updateEkK(oS, k):\n",
    "    #在alpha更新后存储计算得到的误差\n",
    "    Ek = calcEkK(oS, k)\n",
    "    oS.eCache[k] = [1,Ek]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def innerLK(i, oS):\n",
    "    #计算误差\n",
    "    Ei = calcEkK(oS, i)\n",
    "    #找出不满足KKT条件的alpha\n",
    "    if ((oS.labelMat[i, 0]*Ei < -oS.tol) and (oS.alphas[i, 0] < oS.C)) or ((oS.labelMat[i, 0]*Ei > oS.tol) and (oS.alphas[i, 0] > 0)):\n",
    "        #选择j\n",
    "        j,Ej = selectJK(i, oS, Ei)\n",
    "        #存储旧的值\n",
    "        alphaIold = oS.alphas[i, 0].copy(); alphaJold = oS.alphas[j, 0].copy();\n",
    "        #两种情况求边界值\n",
    "        if (oS.labelMat[i, 0] != oS.labelMat[j, 0]):\n",
    "            L = max(0, oS.alphas[j, 0] - oS.alphas[i, 0])\n",
    "            H = min(oS.C, oS.C + oS.alphas[j, 0] - oS.alphas[i, 0])\n",
    "        else:\n",
    "            L = max(0, oS.alphas[j, 0] + oS.alphas[i, 0] - oS.C)\n",
    "            H = min(oS.C, oS.alphas[j, 0] + oS.alphas[i, 0])\n",
    "        if L==H: return 0\n",
    "        #计算变化量\n",
    "        eta = 2.0 * np.dot(oS.X[i:i+1,:], oS.X[j:j+1,:].T) - np.dot(oS.X[i:i+1,:], oS.X[i:i+1,:].T) - np.dot(oS.X[j:j+1,:], oS.X[j:j+1,:].T)\n",
    "        if eta >= 0: return 0\n",
    "        #更新alpha\n",
    "        oS.alphas[j, 0] -= oS.labelMat[j, 0]*(Ei - Ej)/eta\n",
    "        #约束alpha\n",
    "        oS.alphas[j, 0] = clipAlpha(oS.alphas[j, 0],H,L)\n",
    "        updateEkK(oS, j)\n",
    "        if (abs(oS.alphas[j, 0] - alphaJold) < 0.00001): return 0\n",
    "        oS.alphas[i, 0] += oS.labelMat[j, 0]*oS.labelMat[i, 0]*(alphaJold - oS.alphas[j, 0])\n",
    "        updateEkK(oS, i)\n",
    "        b1 = oS.b - Ei- oS.labelMat[i, 0]*(oS.alphas[i, 0]-alphaIold)*np.dot(oS.X[i:i+1,:], oS.X[i:i+1,:].T) - oS.labelMat[j, 0]*(oS.alphas[j, 0]-alphaJold)*np.dot(oS.X[i:i+1,:], oS.X[j:j+1,:].T)\n",
    "        b2 = oS.b - Ej- oS.labelMat[i, 0]*(oS.alphas[i, 0]-alphaIold)*np.dot(oS.X[i:i+1,:], oS.X[j:j+1,:].T) - oS.labelMat[j, 0]*(oS.alphas[j, 0]-alphaJold)*np.dot(oS.X[j:j+1,:], oS.X[j:j+1,:].T)\n",
    "        if (0 < oS.alphas[i, 0]) and (oS.C > oS.alphas[i, 0]): oS.b = b1\n",
    "        elif (0 < oS.alphas[j, 0]) and (oS.C > oS.alphas[j, 0]): oS.b = b2\n",
    "        else: oS.b = (b1 + b2)/2.0\n",
    "        return 1\n",
    "    else: return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def smoPK(dataMatIn, classLabels, C, toler, maxIter):\n",
    "    #建立类变量\n",
    "    oS = optStructK(np.array(dataMatIn),np.array(classLabels).reshape(-1, 1),C,toler)\n",
    "    iter = 0\n",
    "    entireSet = True; alphaPairsChanged = 0\n",
    "    #执行循环\n",
    "    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):\n",
    "        alphaPairsChanged = 0\n",
    "        if entireSet:\n",
    "            #遍历所有\n",
    "            for i in range(oS.m):        \n",
    "                alphaPairsChanged += innerLK(i,oS)\n",
    "                #print(\"fullSet, iter: {} i:{}, pairs changed {}\".format(iter,i,alphaPairsChanged))\n",
    "            iter += 1\n",
    "        else:\n",
    "            #遍历非边界值\n",
    "            nonBoundIs = np.nonzero((oS.alphas > 0) * (oS.alphas < C))[0]\n",
    "            for i in nonBoundIs:\n",
    "                alphaPairsChanged += innerLK(i,oS)\n",
    "                #print(\"non-bound, iter: {} i:{}, pairs changed {}\".format(iter,i,alphaPairsChanged))\n",
    "            iter += 1\n",
    "        if entireSet: entireSet = False \n",
    "        elif (alphaPairsChanged == 0): entireSet = True  \n",
    "        #print(\"iteration number: {}\".format(iter))\n",
    "    return oS.b,oS.alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b = [[-2.89901748]]\n",
      "alphas[alphas > 0] = [ 0.06961952  0.0169055   0.0169055   0.0272699   0.04522972  0.0272699\n",
      "  0.0243898   0.06140181  0.06140181]\n"
     ]
    }
   ],
   "source": [
    "dataArr, labelArr = loadDataSet(\"testSet.txt\")\n",
    "b, alphas = smoPK(dataArr, labelArr, 0.6, 0.001, 40)\n",
    "print(\"b = {}\\nalphas[alphas > 0] = {}\".format(b, alphas[alphas > 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvSSMhhQ4CARI6IQmJIIKIFFFZQLChoCK2\ndVFQ1J/uqiuoa2PXgq6KioqgAiIoioJYEEQElyamUaQECDUE0nvm/f0xEFrKpEzuZPJ+nmeeZMq9\n571DeOfMOee+14gISiml3J+H1QEopZSqGZrwlVKqjtCEr5RSdYQmfKWUqiM04SulVB2hCV8ppeoI\nTfhKKVVHaMJXSqk6QhO+UkrVEV5WB3Cmpk2bSkhIiNVhKKVUrbJp06ZjItKsvNe5VMIPCQlh48aN\nVoehlFK1ijFmryOv0yEdpZSqIzThK6VUHaEJXyml6giXGsMvSUFBAUlJSeTm5lodinKAr68vwcHB\neHt7Wx2KUuocLp/wk5KSCAwMJCQkBGOM1eGoMogIKSkpJCUlERoaanU4SqlzuPyQTm5uLk2aNNFk\nXwsYY2jSpIl+G1PKRbl8wgc02dci+m+llOuqFQlfKaUAiopy+fPPyeTnH7M6lFpJE34N+fLLL0lI\nSCi+P3XqVH788ccaj2PJkiVMmzatzNckJiYyb968GopIKceI2Ni27TYOHHiD9PR1VodTK2nCryHn\nJvx//etfDBkypMbjGDlyJI899liZr9GEr1zR7t2Pk5y8kA4dXqJp06utDqdWcquEHxQExpx/Cwqq\n/D4TExPp2rUrt9xyC926deOGG24gOzsbgBUrVhAdHU1ERAR33nkneXl5ADz22GOEhYURGRnJI488\nwtq1a1myZAmPPvooUVFR7Nq1i9tvv51FixaxfPlyRo8eXdzeqlWrGDFiBADff/89ffv25cILL2T0\n6NFkZmaeF9/AgQOZPHkyUVFRhIeHs379egCOHz/ONddcQ2RkJH369CEmJgaA2bNnM2nSJABuv/12\nHnjgAS655BLat2/PokWLiuP/5ZdfiIqKYvr06ZV/85SqJgcOvMP+/f+hVav7CA5+2Opwai8RcZlb\nz5495VwJCQnnPVYaKP1WWXv27BFA1qxZIyIid9xxh7z00kuSk5MjwcHBsn37dhERGTdunEyfPl2O\nHTsmnTt3FpvNJiIiJ06cEBGR8ePHy8KFC4v3e+p+QUGBtGnTRjIzM0VEZMKECfLxxx9LcnKy9O/f\nv/jxadOmyTPPPHNefAMGDJC7775bRER+/vln6d69u4iITJo0SZ5++mkREVmxYoX06NFDREQ+/PBD\nmThxYnEMN9xwgxQVFUl8fLx06NBBRERWrlwpw4cPr/R7VpF/M6XKc+zYUlm50kNiYkZIUVGB1eG4\nJGCjOJBj3aqH7yxt2rShX79+ANx6662sWbOG7du3ExoaSufOnQEYP348q1evpkGDBvj6+nLXXXfx\nxRdfUL9+/TL37eXlxdChQ/n6668pLCxk6dKljBo1it9++42EhAT69etHVFQUc+bMYe/ekusjjR07\nFoDLLruM9PR0UlNTWbNmDePGjQNg8ODBpKSkkJ6eft6211xzDR4eHoSFhXHkyJFKv0dKOUNGxmbi\n428kICCKbt3m4+Hh8qcOuTR99xxw7lLDspYeenl5sX79elasWMGiRYt48803+emnn8rc/5gxY3jz\nzTdp3LgxvXr1IjAwEBHhiiuuYP78+dUa37nq1atX/Lu9o6CUa8jN3Uds7Ai8vZsQEfENXl4BVodU\n62kP3wH79u1j3Tr7qoB58+Zx6aWX0qVLFxITE9m5cycAH3/8MQMGDCAzM5O0tDSGDRvG9OnT+eOP\nPwAIDAwkIyOjxP0PGDCAzZs389577zFmzBgA+vTpw6+//lq8/6ysLHbs2FHi9gsWLABgzZo1NGjQ\ngAYNGtC/f3/mzp0L2OcFmjZtSpCDkxllxapUTSgsTCM2djhFRdlERi6jXr2WVofkFjThO6BLly68\n9dZbdOvWjRMnTnDvvffi6+vLhx9+yOjRo4mIiMDDw4MJEyaQkZHBiBEjiIyM5NJLL+XVV18F7L34\nl156iejoaHbt2nXW/j09PRkxYgTffvtt8YRts2bNmD17NmPHjiUyMpK+ffuybdu2EuPz9fUlOjqa\nCRMm8MEHHwDw9NNPs2nTJiIjI3nssceYM2eOw8cbGRmJp6cnPXr00ElbVeNstnzi4q4nO3sb4eFf\n4O/f3eqQ3IZxpa/xvXr1knMvgLJ161a6devm0PZBQVBSxzQwEEoYvnZIYmIiI0aMIC4urnI7cLKB\nAwfy8ssv06tXL6tDKVaRfzOlziQibN9+J4cPz6Zr19lccMF4q0OqFYwxm0Sk3CTgVmP4lU3qSinX\nsHfvcxw+PJuQkKc12TuBWyV8ZwgJCXHZ3j3Yx+eVcgeHD39MYuJUWrQYT7t2U60Oxy3pGL5SynIn\nTqxk+/a7aNhwMF26zNQifE6iCV8pZamsrATi4q7Fz68z3bt/joeHj9UhuS1N+Eopy+TlHSYmZhie\nnn5ERi7F27uh1SG5NR3DV0pZoqgoi9jYERQUJBMdvRpf33ZWh+T2tIfvAlJTU5kxY0bx/YMHD3LD\nDTdYEsuwYcNITU0t8zWzZ8/m4MGDNRSRckciRSQkjCUz83fCwhYQGNjT6pDqBE34LuDchN+qVavi\nypU1bdmyZTRsWPbXak34qipEhJ07HyQl5Ws6dXqDpk1HWB1SnaEJvxyJiYl069aNv/71r3Tv3p0r\nr7ySnJwcAHbt2sXQoUPp2bMn/fv3Lz4TdteuXfTp04eIiAiefPJJAgLsNUAyMzO5/PLLufDCC4mI\niOCrr74C7OWId+3aRVRUFI8++iiJiYmEh4cD9hIL8fHxxfEMHDiQjRs3kpWVxZ133knv3r2Jjo4u\n3teZVq1axWWXXcbw4cPp0qULEyZMwGazATB//nwiIiIIDw/nH//4R/E2ISEhHDt2rNTjXrRoERs3\nbuSWW24hKiqq+L1QylFJSa9x4MCbtGnzCK1b32d1OHWLIyU1a+pW1fLIxdauFXnhBfvPKtqzZ494\nenrK77//LiIio0ePlo8//lhERAYPHiw7duwQEZHffvtNBg0aJCIiw4cPl3nz5omIyNtvvy3+/v4i\nIlJQUCBpaWkiIpKcnCwdOnQQm80me/bsKS5rfKrNU/dfffVVmTp1qoiIHDx4UDp37iwiIo8//nhx\nHCdOnJBOnToVl1I+ZeXKlVKvXj3ZtWuXFBYWypAhQ2ThwoVy4MABadOmjRw9elQKCgpk0KBBsnjx\nYhERadeunSQnJ5d53AMGDJANGzaU+p5peWRVmqNHF8nKlUbi4m4Qm63I6nDcBg6WR3a/Sdt16+Dy\nyyE/H3x8YMUK6Nu3SrsMDQ0lKioKgJ49e5KYmEhmZiZr16496+Ilpy6Asm7dOr788ksAbr75Zh55\n5BHA/uH6xBNPsHr1ajw8PDhw4EC5JYlvvPFGrrzySp555hk+++yz4rH977//niVLlvDyyy8DkJub\ny759+84radC7d2/at28P2Msor1mzBm9vbwYOHEizZs0AuOWWW1i9ejXXXHNNucetVGWlpf3G1q23\nEhTUh65dP8IYHWCoae6X8Fetsif7oiL7z1Wrqpzwzywh7OnpSU5ODjabjYYNG7JlyxaH9zN37lyS\nk5PZtGkT3t7ehISEkJubW+Y2rVu3pkmTJsTExLBgwQLeeecdwP7h8fnnn9OlS5cyt6+u0smnjlup\nysjJ2UVc3NX4+LQmPPwrPD39rA6pTnK/j9iBA+09e09P+8+BA53STFBQEKGhoSxcuBCwJ+BTpZD7\n9OnD559/DsCnn35avE1aWhrNmzfH29ublStXFl/QpLxyxDfddBP/+c9/SEtLIzIyEoCrrrqKN954\no7iG/e+//17ituvXr2fPnj3YbDYWLFjApZdeSu/evfn55585duwYRUVFzJ8/nwEDBjh87Fo+WVVE\nQUEKMTF/QUSIjPwWH59mVodUZ7lfwu/b1z6M8+yz1TKcU5a5c+fywQcf0KNHD7p37148cfraa6/x\n6quvEhkZyc6dO2nQoAFgHzrZuHEjERERfPTRR3Tt2hWAJk2a0K9fP8LDw3n00UfPa+eGG27g008/\n5cYbbyx+bMqUKRQUFBAZGUn37t2ZMmVKiTFedNFFTJo0iW7duhEaGsq1115Ly5YtmTZtGoMGDaJH\njx707NmTUaNGOXzct99+OxMmTNBJW1WuoqJc4uKuITd3HxERX1G/fierQ6rT3Ko8sqvIzs7Gz88P\nYwyffvop8+fPL3EVjbOtWrWKl19+mW+++aZG262N/2aq+onYSEi4meTkBYSFLaB58xvL30hVSp0s\nj+wqNm3axKRJkxARGjZsyKxZs6wOSakat2fPP0lOXkD79v/WZO8iNOE7Qf/+/YvH8600cOBABjpp\nDkOpshw8OJN9+6bRqtUE2rQ5f5hSWcP9xvCVUpZKSfmWHTvuo3HjYXTs+IaWOnYh1dLDN8bMAkYA\nR0Uk/ORjjYEFQAiQCNwoIieqo73SbD60mXc2vsMfR/4gpyCHRn6NGN5pOHdG30nT+k2d2bRSCsjI\n2EJCwo0EBEQSFrYADw8dRHAl1dXDnw0MPeexx4AVItIJWHHyvlNsO7aNvh/05doF1xLaMJTXh77O\nx9d+zJP9n2TbsW10eqMT9y29j7zCPGeFoFSdl5u7n9jY4Xh5NSIi4hu8vAKsDkmdo1o+fkVktTEm\n5JyHRwEDT/4+B1gF/INqtuXwFq765CqeHfQsd0XfhaeH51nPX9HhCl6+8mXu+foehs0bxtKbl+Lr\n5VvdYShVpxUWphEbO5yiokyio9dQr14rq0NSJXDmGH4LETl08vfDQIvqbiA1N5Wr51/NW8Pe4p6e\n95yX7E9p7NeYBTcsoLFfY+5fdn91h1FjXnvtNbKzs522/y1btrBs2TKn7V+5J5utgPj40WRnb6V7\n988JCIiwOiRVihqZtD1Z3KfEBf/GmHuMMRuNMRuTk5MrtN85W+bQv21/bggrv3a8p4cnH4z8gM+3\nfs7BjNpX2reoqMglE35hYaGTolG1gYiwY8cETpz4gc6dZ9K48RCrQ1JlcGbCP2KMaQlw8ufRkl4k\nIjNFpJeI9DpVzMsRIsKMjTOYeNFEh7cJqhfEmPAxvLfpPYe3ycrKYvjw4fTo0YPw8HAWLFgAnC4j\nDLBx48bi5Y9PP/0048aNo2/fvnTq1In33rO3VZlSxQEBAfzf//0fPXr04Pnnn+fgwYMMGjSIQYMG\nnRXj8uXLzyritmrVKkaMsNcY//777+nbty8XXngho0ePJjMzE4ANGzZwySWX0KNHD3r37k1aWhpT\np05lwYIFREVFsWDBAo4fP84111xDZGQkffr0ISYm5qxj7NevH+PGjXP4vVTuZ9++Fzh8eBbt2k2h\nZcs7rA5HlceRkpqO3LCvxok74/5LwGMnf38M+E95+6hIeeStyVsl5LUQsdlsjtUPPenXfb9K9DvR\nDr9+0aJFcvfddxffT01NFZHTZYRFRDZs2CADBgwQEZGnnnpKIiMjJTs7W5KTkyU4OFgOHDhQqVLF\ngCxYsKC47TPbPFNBQYG0adOmuDzyhAkT5OOPP5bk5GTp379/8ePTpk2TZ555RvLy8iQ0NFTWr18v\nIiJpaWlSUFAgH374oUycOLF4v5MmTZKnn35aRERWrFghPXr0KD7GCy+8ULKzs0t8z7Q8ct1w+PAn\nsnIlkpAwrsL/D1X1wsHyyNXSwzfGzAfWAV2MMUnGmLuAacAVxpg/gSEn71eblOwUWga0rPAa35YB\nLTmec9zh10dERPDDDz/wj3/8g19++aW4Lk5ZRo0ahZ+fH02bNmXQoEGsX78eOF2q2NPTs7hU8YYN\nG4pLFXt5eRWXKgZ7hcrrr7++3Pa8vLwYOnQoX3/9NYWFhSxdupRRo0bx22+/kZCQQL9+/YiKimLO\nnDns3buX7du307JlSy666CLAXgjOy+v8+fs1a9YU9+AHDx5MSkoK6enpAIwcORI/P614WFelpv7M\ntm130LDhQLp0eV/X2tcS1bVKZ2wpT11eHfsvia+XL7mFZZcWLkluYW6FVul07tyZzZs3s2zZMp58\n8kkuv/xypk6dipeXV/GQzLkljksrSVzRUsW+vr54epY8EX2uMWPG8Oabb9K4cWN69epFYGAgIsIV\nV1zB/Pnzz3ptbGysQ/ssi7+/f5X3oWqnrKytxMVdg59fR7p3/wIPDx+rQ1IOqrVn2rZv1J49qXtI\nyU6p0HZr9q2hWzPHC3sdPHiQ+vXrc+utt/Loo4+yefNmwD6Gv2nTJoDiUsinfPXVV+Tm5pKSksKq\nVauKe9JVLVVcVlniAQMGsHnzZt577z3GjBkD2Ms0//rrr+zcuROwz0fs2LGDLl26cOjQITZs2ABA\nRkYGhYWF5+2/f//+zJ07F7DPCzRt2pSgoCCH3zvlfvLzjxAbOwxj6hERsQxv70ZWh6QqoNYm/EZ+\njRjZZSQfbvnQ4W3k5ETvhJ4THN4mNjaW3r17ExUVxTPPPMOTTz4JwFNPPcXkyZPp1avXeb3wyMhI\nBg0aRJ8+fZgyZQqtWtnXJFe1VPE999zD0KFDz5u0Bfvwz4gRI/j222+LJ2ybNWvG7NmzGTt2LJGR\nkfTt25dt27bh4+PDggULuP/+++nRowdXXHEFubm5DBo0iISEhOJJ26effppNmzYRGRnJY489xpw5\ncxx+35T7KSrKIjb2avLzjxIR8Q1+fiFWh6QqypGB/pq6VfSatr/t/03aTm8rKdkpDk1sLNm2RDr+\nt6MUOfFamk899ZS89NJL5z2+cuVKGT58uNPadSU6aet+bLZCiYkZJStXekhy8ldWh6POQU1O2lrl\n4uCLGR02mhHzRnAip+wyPav3rubOJXcye9RsPPRamkpVyM6dD5OS8hUdO75O06YjrQ5HVVKtvwCK\nTWw8+v2jLN62mIf7Psy4yHE08D29kmbbsW28veFt5sfNZ9718xjSXk8McTa9AIp7SUp6nZ07HyQ4\n+GE6dnzF6nBUCdzqAigiUuqKFg/jwStXvcLILiOZsXEGU1ZOIbJFJH5efhzNOsrBjIPcfeHdbP7b\nZoKDgms48rrHlToQquqSkxezc+dDNG16HR06vGR1OKqKXD7h+/r6kpKSQpMmTcpcxjggZAADQgZw\nOPMw245tI7sgm0a+jejZqic+nrpsrCaICCkpKfj6anE6d5Ce/j+2br2FoKCL6dbtE4wOhdZ6Lp/w\ng4ODSUpKoiJ1dlqcqtOWCbt27HJSZKokvr6+BAfrN6naLidnN7GxV+Pj05Lw8K/w9NST7NyByyd8\nb29vQkNDrQ5DqTqjoOA4MTHDECkiMvJbfHyaWx2SqiYun/CVUjXHZssjLu5acnP30KPHj9Sv39nq\nkFQ10oSvlAJAxMa2bXeQlraabt3m07Bhf6tDUtVMZ2GUUgDs2TOFo0fnExr6Ii1ajLE6HOUEmvCV\nUhw8+D779r1Ay5b30LZttV+JVLkITfhK1XHHj3/Hjh0TaNx4KJ06vaWljt2YJnyl6rDMzD+Ijx+N\nv384YWGf4eGh03ruTBO+UnVUbm4SMTHD8fJqQGTkUry8Aq0OSTmZfpwrVQcVFqYTGzucoqJ0oqPX\nUK9ea6tDUjVAE75SdYzNVkB8/I1kZcUTGbmMgIBIq0NSNUQTvlJ1iIjw55/3ceLEd3Tp8j6NG19p\ndUiqBukYvlJ1yL590zh06H3atv0nLVveZXU4qoZpwleqjjhyZD579jxB8+Y3Exr6rNXhKAtowleq\nDkhNXc22bbfToMEAunadpWvt6yhN+Eq5uezs7cTFXYOfX3vCwxfj4VHP6pCURTThK+XG8vOPEhPz\nF4zxJiJiGd7ejawOSVlIV+ko5aaKirKJjR1Jfv5hoqJW4een15Wo6zThK+WGRIrYuvVWMjLWEx6+\nmKCg3laHpFyAJnyl3NCuXY9y7NhiOnZ8naZNR1kdjnIROoavlJtJSnqDpKTptG49meDgB6wOR7kQ\nTfhKuZFjx75i587JNG16DR07vmJ1OMrFaMJXyk2kp28gIWEsgYEX0a3bXIzxtDok5WI04SvlBnJy\n9hAbOwIfnwuIiFiCp2d9q0NSLkgTvlK1XEHBCWJjhyNSQETEMnx8WlgdknJRukpHqVrMZssjPv46\ncnJ20aPHD/j7d7U6JOXCNOErVUuJCNu3301q6iq6dZtLw4aXWR2ScnE6pKNULZWY+BRHjnxCaOjz\ntGhxs9XhqFrA6T18Y0wikAEUAYUi0svZbSrl7g4dmsXevc/SsuXdtG37uNXhqFqipoZ0BonIsRpq\nSym3dvz4D+zY8TcaNbqSTp1maKlj5TAd0lGqFsnMjCU+/nrq1w+je/eFeHh4Wx2SqkVqIuEL8KMx\nZpMx5p5znzTG3GOM2WiM2ZicnFwD4ShVO+XlHSA2dhienoFERCzFyyvI6pBULVMTCf9SEYkC/gJM\nNMactZRARGaKSC8R6dWsWbMaCEep2qewMIPY2BEUFqYSEbEUX99gq0NStZDTE76IHDj58yiwGNA6\nrUpVgM1WSELCjWRmxhIWtpDAwCirQ1K1lFMTvjHG3xgTeOp34EogzpltKuVORIQ//5zI8ePL6dz5\nbZo0GWp1SKoWc/YqnRbA4pOrCLyAeSKy3MltKuU29u//D4cOzaRt28dp1eqvVoejajmnJnwR2Q30\ncGYbStWUzPxMlu9czpHMI3h7etOhUQcGhQ7Cwzjni/LRowvYvfsxmjcfQ2joc05pQ9UtWlpBqXLs\nPrGb6eumMzd2Ln2C+9CuQTsKbYW8s/Ed0vLSmNBzAvdddB/+Pv7V1mZq6hq2br2NBg3607XrbIyT\nPlRU3aIJX9U+69bBqlUwcCD07evUplbvXc0Nn93AXy/8KzH3xhAcdHp1jIiw/sB6/v3rv/k0/lOW\n3byMFgFVr1SZnb2DuLhR+PqGEB6+GA+PelXep1IARkSsjqFYr169ZOPGjVaHoVzZunVw+eWQnw8+\nPrBihdOS/h+H/+CKj69g3vXzGNJ+SKmvExGe/OlJlu9azurbV1epp5+fn8zmzX0oKsrgwgt/w8+v\nfaX3peoOY8wmR8rW6PdEVbusWmVP9kVF9p+rVjmtqYe/f5jnBj9XZrIHMMbw3ODnaNugLe9sfKfS\n7RUV5RAXN5L8/IOEhy/RZK+qnSZ8VbsMHGjv2Xt62n8OHOiUZrYmbyX+aDzje4x36PXGGP5+yd95\ne+Pb2MRW4fZEbGzdOo709P/RrdtcGjToU+F9KFUeTfiqdunb1z6M8+yzTh3OmfX7LO6MvpN6XiWM\nn69bBy++aP95hj7BfQjwCeDnxJ8r3N6uXX/n2LHP6dDhFZo1u66yYStVJp20VbVP375On6zddWIX\nN0eUUGO+jDkEYwzRLaPZfWI3g0IHOdzWgQNvkZT0Cq1b309w8IPVdQhKnUd7+EqVoNBWiJdHCf2h\ncuYQvD28KbQVOtzOsWNf8+efD9CkyUg6dpyupY6VU2nCV6oELfxbkJiaeP4T5cwh7EndQ3P/5g61\nkZ6+kYSEMQQGXkhY2DyM8axy3KpygoLAmPNvQW5WkFQTvlIlGBM+hlm/z+K8ZctlzCEkpiay+dBm\nrup4Vbn7z83dS2zsCHx8mhMR8Q2entV30paquIyMij1eW+kYvlIlGBw6mPyifFbvXc2AkAFnP1nK\nHMI7G9/htsjbqO9dv8x9FxSkEhMzDJE8IiJW4uNT9ZO1VNmCgkpO3oGBkJ5e8/FYRXv4SpXAGMMz\nA5/h9q9u52DGwXJfv3zncmZvmc3kPpPLfJ3Nlk98/HXk5PxJ9+6L8ffvVl0hqzLUlR58eTThq1ol\nLTeNncd3sufEHrILsp3a1k3hN3Ffr/voN6sfqxJXnT+8A+QW5jJjwwxuW3wbn9/4OSENQ0rdn4iw\nffvdpKaupEuXWTRqNNB5wStVAh3SUdaoQD0cm9j4ftf3zNgwg5WJK2nu35wiWxEnck9wfbfrue+i\n++jVqtyzyivl0X6P0q5hO+5dei+expM7ou4gpGEIBbYCNh/azOwts4m6IIqfxv9EePPwMveVmPgM\nR458TEjIs1xwwa1OiVepsmgtHVXzKlAP50TOCa5dcC1peWlMumgSYyPGFo+RJ2cl8+GWD5mxYQZD\n2g/h7eFv4+3pnIt6iwirElexIH4BR7OO4u3pTfuG7bkz+k46NelU7vaHDs1m+/Y7uOCCO+nS5X1d\nflnDynq7RWr/GL+jtXQ04aua9+KLMGWKfS27p6d9xcvjj5/3ssz8TAbOHkj/tv155apXSq07n5mf\nyY0Lb6SBbwPmXjfXafXpK+vEiRXExAylYcOBREQsw8PDOR9KqnTlJfzaTounKdflYD2cf674J2HN\nwnj1qlfLTOIBPgF8fuPn7D6xmzlb5jgn5krKzIwjLu466tfvRvfuizTZVzMRYd3+dczZMof3Nr3H\n4q2LyczPPO91gYElb1/a4+5Ke/jKGuWM4WfmZ9J2etvzatCXtY/lO5fzxIon2HTPJpcYMsnLO8jm\nzX0QKeTCC/+Hr28bq0NyG3mFecz6fRYzNs6goKiA3q174+Ppw/70/Ww4sIFbIm7hwT4P0qFxh7O2\nq66hG1cbAnK0h6+Ttsoa5dTDmRc7j4EhA0tO9uvWwUcfwYcfQmFh8TzAlX2uZOKyiaw/sJ6Lgy92\nYvDlKyzMJDb2agoLTxAVtVqTfTVKzU3lmk+vwcfTh5eveJlL2lxCgE9A8Yf8/rT9vLvpXfp+0Jf5\n18/n8vaXF29bXcsza+syTx3SUS5pw4ENXNWhhDNWT034vvsu5OWdVdPGw3hwVYer2HjQ2m+JNlsh\nCQljyMz8g7CwzwgMjLY0HneSX5TP8HnDySnIYfux7YyYP4IWL7eg/gv1Gbd4HOv2ryM4KJjnBj/H\nwtELGfv5WMv/HlyJJnzlkrIKsgjwCTj/iVPFy04NRRpz1jyAv7c/WQVZNRbnuUSEnTsf4PjxpXTu\n/BZNmvzFsljc0fjF41l/YD0RLSL4csyX5D+ZT/Y/s0l6KInoC6K57cvbuOqTq0jNTWVAyABevepV\nJi2bZHW2L30fAAAgAElEQVTYLkOHdFTVOOn6skH1gkjJSTn/iVMTvvn59knfYcPggguKnz6ec/y8\ncduatH//Kxw8+DZt2vyDVq3+Zlkc7ujFX15k8bbFvH/1+4yPsl+Y5vRYehPgYTCT2Tn0IRqtuRxm\n/wwFY2HyPzGtNsOhC60M3yXopK2qmDMTPDjt+rKfJ3zO6/97ndV3rC49hiZN4MEHi9vP/X4Zbdfd\nyK93/urQ2vjqdvToQhISbqRZs5tOVr/UL9DVZc2+NVz/2fX4evmyZ/Ke4lVbJc/NC4y8GzCw5H3o\n/wIEJcHSGWW2UZFUWNVlntU96auTtqriyuutn3vC1Pjx59eGdzDhF9oKWbpjKeuS1pGZn0lQvSAG\nhQzi8vaX42E8GNllJA8sf4DYI7FEtIg4e+NTE74vvnhW+4tWzSA6JNqSZJ+WtpatW8cRFNSPrl1n\na7KvZv/9338Z0WkEaXlpDpxnYeCH/8ADHeHHaXCwJ4SsKnOLii7PDAwsfYLWmPITd1mTvo5sX1n6\nV6nsTiXzKVPsP8+5fB9gXxmTm3s6wUOFry9baCtk2ppptH+9Pf/+9d8E+gTSuUlnfL18eeSHR+j6\nZldmLHgEr3+/xAMXjGLisonkFeaVvLMz1vMfbuTNFO81PNTnoUq/BZWVnf0nsbEj8fVtS0TEV3h6\n+tZ4DO7sUMYhftj9A5e0uaT4ojSn6teXKqcJbB8F0bPA5g2e+SW+TMR+q2hyTU8vuydf1dU6zlrt\noz18ZVfSlZzO7K2vWwezZp3+K/fygttus98cHMPPLczlugXXUWgr5KsxXxHd8uzVK1Mum8K65e/x\n0JcT+O0YfPC8D5te6cfV86/ms9Gf0dC34dk7PFmbfvfKLxjus5C7et/N0I5Dq/pOVEh+/jFiY4dh\njCEiYhne3k1qtP264Lek37i07aW0a9iOfWn7AAcT4rZREP0hZDWDLMcuSuPuNOEruzMnQ0vqra9a\nZf8wAHvX6o47Tid4B4dx7vjqDgJ8Aph73dwSa94YY7hkSwor5xhGjLHxj8vymXd8EA+3O0r719tz\nc8TN3H3h3bRv1J4iWxExR2J45+A7LPdczr/6/4v7L76/0odfGUVFOcTFjSIvL4kePX6ifv2ONdr+\nuU4Nk/24+0fS8tLw9/and+ve3BR+U7k1+l1ZRn4GDX0b0r9tf3af2M3W5K2AA2Wl8xpAvXR7L/83\nvVYw6JCOOqWMKzkBZ5dD8PW19+wrYMOBDazbv46Prv2o7AJnAwdS36MeXyz0YE4P4UCfMP77l/8S\nc28MTfyaMHrhaFq/2prQ10N5YPkD9A3uS+LkxBpP9iI2tm0bT3r6Orp1+4QGDZx7UfWy2MTGK2tf\nKR4mC20UyuDQwXRv3p0vtn1B2+ltefT7R51eTtpZAn0CSctNo1njehz59i7Cxpc9+VqsXhog0Gg3\nbB9Z4Xbd8rKHIuIyt549e4pyYWvXirzwgv1nRZ4Tkdu/vF2m/TKtQu1MnnWjPPHjE1UI2Hl27nxU\nVq5E9u172dI4CooKZMyiMdLvg36y+eDmEl+z+/huGbNojPR5v4+k5qTWcIRVdyD9gDSa1kiolyYE\nHhD+r6XQ5Us5PQJf8s3zprHS/KXmUq/f2yU+HxhYdrtl7bsirylJYGDZ25a3/fmxslEcyLGWJ/kz\nb5rwa5EzE/zatSJ+fiKenvaf5yT9/MJ88X/eX45mHi17P+eIPxovwa8GO+sIKi0paYasXIls3z5R\nbDabpbFM/nayDPloiOQU5JT5OpvNJvd+c68M+WiIFNmKaii66nPDZzcIF71lT4at1guPNBd6/1fw\nyik5YYb+IEz1kAe+faDSbTqSjEtL3OV9mFTX9qdjdSzh6xh+XVDOcsuU7BQ+3PIhs7fMJik9CUFo\nFdiKcZHjuCv6LloEtDh/f6eWZ3p6QlSUvcyBzVbihO/xnOPU965PM/9mpe+nhHX8nZt05lDGIWxi\nc5mSxykpS/nzz0k0aTKCjh1fs7RIW1J6Eh/98RG7J+/G16vslUHGGN74yxtEvRvFT3t+Ykj7ITUU\nZfV4oPcDLLr0Fth6LRy8yH5S1dDJMOBZ+P0OONAbinzs6+2jZkOzBEjqy+tDX3dqXFVdOlnThdY0\n4bu70pLqunXIypVMa7eff++fz6iuo5h59UzCmoUBsCNlB+9teo+ub3Xlvl738ezgZ+1Jd906ePrp\n0wm+qAjWr7e35eFR6vJMoYQ1bOWsDLJ3XFxHRsZm4uNvIiAgmrCwT/HwsPa/z8xNM7kl4pazVi+V\nfUKPJxMvmlh8wRjHtql6nJsObmLGhhn8sPuHsyaT7+11L1d0uKLUD/Oz4+oP/SfAbUNg7rdwrCt8\n8h00/hMufB8i54JnHmQ3gcwL7L9/8i1BC2vHBUxqjCNfA2rqpkM6TvDCC/ahFrD/PDV84ucnj15p\nJOpeI/t/+rLUzY9kHpFLPrhE/rrkr2L79Vf7kI2Hx/nfQT08RK68ssShmfzCfAl8IVAOZRw6+4ly\nhoJiDsdI2+ltq+VtqKqcnL3y668XyNq1bSU396DV4YiISPvX28uWQ1vOeqy8YYj03HTxf95f0nPT\nHd7mTMlZyfLvNf+W3u/1lvavt5eub3aVUfNHydIdS88aKgoMFCHgkHD7AOHBdsKlLwqN/5SApsdl\nf9p+qXfJe8LfooRJnYULfi9xKKPEmPq8KvyjoTDinpPb2eyP+x4X+ky37++2wfb7pRyDo8p6Xyo6\n5OJs6Bi+krVrRSZMEKlX7+yk+sILsiDCSJdJSIq/h/1DoAwZeRkS/UpnmXlr2Olk7+Eh0rv3+fsu\nxd1f3S3P/fxc6TFOmHDe9hOXTpSpP02t1KFXp4KCVPnf/7rL6tUNJDMzzupwitV/vr5k5GWc9Zgj\nybvd9Hay6/iuCm1TUFQgDy1/SBpOayi3f3m7rNyzUnYc2yHxR+Plg80fSM93e0roa6Hy0+6f7PsM\nOCg80F4Y8LRgCs/bp/13mxA+X3i0mdBqg+MTogEHhcv+JTwULEzxFp7wF56sJ1x3s9BmzekPgSom\n/PImVl2JyyR8YCiwHdgJPFbWazXhl6GcVTAlvv5U79nH5+yEunat9L7HyNddPcpN1Kdev7qTj3SZ\nhNhOJftT2zkY1+aDm6XNq20kKz+r9DjPiOVY1jFpNK2R7E/b79jxOklRUZ78/vvlsmqVtxw/vsLS\nWM4V8EKApOWmnfWYIwkq+NVgSTyR6PA2BUUFcs2n18jQT4ZKclZyibHYbDZZ/udyafafZrJk2xLh\njkuFAc+Uus+zHuvypfBwK8Enw7GEf+bNK8e+escUOfT6ik+GulfCd+pMmDHGE3gL+AsQBow1xoQ5\ns0235EjZg3OdOT5eVARt2xaPj29s582R9i34y7h/nTWmz4svlrzvVau4dHch3kWwsr2BIUNOb9e3\nr/16tOWcfBXdMprBoYMZ+/lY8ovOOM29hHH8jLwMRn06ir/1/FvpV7uqASLCjh1/IzV1BV26vE+j\nRoMti6Uk7Rq0I/ZIbIW2Sc5KJj0v/fwJ9DJMXTmVrPwsloxZQtP6TUt8jTGGqzpexTc3f8NtX94G\nDRNh9ZOONbB9lH3SNWKewzEVK/SFvCAQx1KZq1+gxNmcvfShN7BTRHaLSD7wKTDKyW26n5ImN8tT\nxnVjf9z9I9dFjcXziX/aH7j3XvvzU6bAoEH2+zNnnv4AGDgQ41OP0VsN33f2tE/aVqIq5syrZ+Jp\nPBny0RDWfjsTeeEFe8XLk3HafLz5IaI+l82+jLBmYTx/+fOnNy7rA8lJ9u59lsOHZxMS8jQXXFCx\nE81qwp3Rd/LOpncqtM2s32dxfbfrHT7zNjM/k3c2vsMHIz8o+4S5k3q37k2boDaQGlJqEi5xYdOG\n++CitwCp/Sc3uTBnLzNoDew/434SYO2151yNI/Xkyyt7UJJTZ86eue+TbaW13kqTkC6nvznk5tq/\npYL9Q+Xdd+33PTygXj37flasoMmKFznc0qvSJZB91m9i4a6evOUXy/gtE6hfIIz42YvA567mRFYK\ni3334Jc0iwcvfpDbo24/veSxnOWbznD48EckJj5FixbjadduqlPbqqzbo27nX6//i31p+2jboC1Q\nehXHwEDIys/i7Y1vs3D0wvOeK22bT2I+YWDIQNo0OH2JxvJW9RzOPAzN94N3FhT4O3Ywey6HBvvB\n/yhktXBqT/zUn5VV15+1kuXLMo0x9wD3ALRt29biaGqYo4mstOT90Uf252+7rfTtTj1+Rlv+AwxZ\n42+GdUlnJ/tTTt0/c13944+TWdQP/+xjVTpWz/x8HvDwYFIRrAiB39oWcTzrGEEDrmRWyED6tel3\n/tr28gq7VbMTJ35i+/a7adhwMF26zHSJC6KXpLFfY6ZcNoUR80aw4rYVNPNvVmoCyy3MZfTCMQwK\nHcRFrS8667mykt7V85cyvsf4sx4r73quGfkZcLQXtF0Du0q4TGVJxMO+pNI3DbJalPvyssoTO8qR\n7cv6MKyNnJ3wDwBnXr05+ORjxURkJjAT7BdAcXI8rqUiiezc5D1okH0tPNirWJa17Tlr53slefDE\njmU8PysNcyq5e3nBiBGwbJn9wuA223nr6r/b9R1/61nJqzideawieHh4csVe4YpDPvDCv8tO4JX5\nhlNJWVkJxMVdh59fZ7p3/xwPDx+ntVUdHu77MKm5qfT9oC/PD36ea7tdi4/n6ZhtYuOnPT8xZeUU\n2gS14d0R71Zo/ydyTtDCv/wEfCZ/b39ycxqD34kKbYdPFuSf/kZQVrI990PKWZ/J7vYNwNkJfwPQ\nyRgTij3RjwFudnKbFeeky/SV68xE5uUF+/bZYykrca9aZX9d/hkTnwUF9t5+Scdwqmd/6kQpDw+u\n3O/DRJPP+haFXHxqwE0E/v53++3U1aRSUoqT69YXHmKrieHaW6+t+rH6+MBrr53e/7nxnnscJX3D\ncYK8vMPExAzD09OPyMileHs3LH8jixljeHbws1zU+iKm/zadB797kJGdR9KkfhMy8zP5btd31POs\nx/297+euC++q8BnLvl6+5BTmVGibnq168v3ufVDg5/hGzWNBzFm9+9KSbbm18J3M2SeqOZUjS3mq\ncgOGATuAXcA/y3qtJcsyyzn5p0banzDBvnSyrBjOXWbp7X16fZiXl309vIeH/fF33z293ZknXp1x\nctQr8+6Xy8cbyffg9HMlrcdfu1aK/HzlmjHI1Cu8qvb+lLeE08J/i8LCTNmwoaf8/LO/pKdvqrF2\nq1v80Xj572//led+fk6mr5sua/auqVK9n4eWP3ReAbvyliouiFsgTPESGu10bGklIgy/Vxj4lEPL\nHh3epwO3ynDFpZq4Si0dEVkGLHN2O5VWw+PD5+nb93St+bJiODNOgL/+9eznZ8609+BtNpg4ESIi\n7Ps4t2d9coXNAxdfxIq9K7n1hng+Wgz1vOrZX3tOD7tw5QomXJlHih888bOtau/PmcNSJbHo30Kk\niISEsWRm/k54+FcEBtbei12HNQsrLo9RHf7W828MmD2AqQOm2v9GHJBdkA3iBc3j4YQDF5RvuhXC\nP4UZp5eY1tYxcpfnyKdCTd1qbQ+/oidFVSaGsl6zdu3ZPf5ze+ulxJedny03vD1IOj3dRKbPnywn\nfv6uuI30Br4yY8Ej0v2lUPnLOA9J93PwJK2qsKCHb7PZZPv2ibJyJZKU9JbT26uNhnw0RF5Z+0rx\n/TLPQPXOFO7rLvT9j5i/N5Xvdn5XZo844WiCtJ3eVuZsmeNwPBXpwQcGVl9FSkfatwqucqZtRW6W\nnWlblYRdXUnKkRjKes2779qHdjxKScylbGuz2eSXvb/I2EVjJfDpehI6GWn/ABL4OHL9M91lxe4V\n9ho6VflAq4iqfnhW0L59r8jKlcjOnY/USHu10a7ju6Tlyy3loy0flfh8ccLzSRfGXSGMukPAJrT9\nRVq81EK4fozQdrWcVfKgWbzwl0nS+N+NZfbvsysUj7OGaqqjfas4mvCN/bWuoVevXrJx40arw6iY\nF1+0n7BUVGQ/yenZZ+1nnlqhtMnnc5d/ljJhmvHLjxy9aQRSUEDzAh+Cvv2pZoe3alhy8ufEx4+m\nWbPrCQtbgHGREsyuKP5oPMPnDad3695MvGgil7W7rHi5qvFNgx4fw8Wvw94B8M3b9guHA8ezT9By\n6EfkRc6wL7nMbgI+meCZj0/c3fy54K/F5xA4ypEJW2emtbLatyqdGmM2iUivcl+nCb+KLDgpqEJO\nLcn88cfTSy09Pe2/lxSvVSuWalha2jr++GMwAQHR9OixAk/PCqwoqaPSctP4OOZjZmyYQXZBNsFB\nweQV5bFx107YfYX9bNnEAcDpjHgqvYgI+9P3k56Xjr+3P8FBwQ6duVvaipjyODOtueIqHU34NclV\nk2QJSzLx8Dg9uWv1NxKLZGfv5Pff++Ll1ZDo6HX4+JRcH0aVTETYdmwbx7KPUc+rHhd3DoXskmvz\nVDW9VHb5pQultRrhaMK3/Exbt1De6hOrnFr1cirZDxkC118PDz5YIycxuaKCghRiY4chIkRELNNk\nXwnGGLo163b6gUpcG92ZvWRd4VM6TfjurJQlmUREuOY3EicrKsolNnYUubn7iIpaQf36nawOyS1U\npvxAeeUZKquu9ewrShN+beXIMFJpZ6ha8Y3E4mEvERvbtt1OevqvhIV9RoMG/Wo8Bnd1bo/8VO89\nI+PsIZlacSaqm3OfhO+q4+jOUJGJYlcYbnKBie3du58gOXkB7dv/h+bNR9do23WNs3rvqurcYx1a\nZS4QUptVpj5+aWqiznx1xlsJBw++y/79/6ZVqwm0afNIjbatao6O3ZfPPXr4VpdHqGnVVT2ypnre\nNVjt8lwpKd+yY8dEGjceRseOb7hsqWNVcTpeX3HukfAtTCiWqK7qkTX1QVlD1S7PlZHxO/HxowkI\niCQsbAEeHu7x5+4O3K3OfG3hHv8DLEoolqqOsfma/KCs4bmE3Nz9xMaOwNu7MRER3+DlFVBjbZfH\nFU/cqaqKniBVW4+ztnOPhA+uMTlZ27jpB2VhYRqxscMoKsokOvpX6tVrZXVIZ3HHSU1HYw8Kqliy\n128C1ct9Er6qHDf7oLTZCoiPH0129jYiIr4lICDc6pDUGSr6oabfBKqXJnzlNkSEHTsmcOLED3Tp\n8iGNGw+xOiSlXIp7LMtUCti793kOH55Fu3ZTadnydqvDUcrlaMJ3BTWxFt7NHT78CYmJU2jRYhwh\nIU9bHY5SLkmHdKzmAmeh1nYnTqxi+/Y7adhwEF26vO/ya+3dcSKytGNSrkV7+Faz+CzU2i4rayvx\n8dfi59eJ7t2/wMPDx+qQypWeXvL1kmrzBOW5x1Tah1dt/lBzB9rDt1pdO2msGuXnHyE2dhjG1CMi\nYine3g2tDkmdVJs/vNyZJnyruelaeGcrKsoiNnYE+flHiYr6GT+/EKtDUsrlacJ3BW62Ft7ZRIpI\nSLiZjIzNhId/SVBQuRf6UUqhCV/VQjt3PkxKyhI6dXqTpk2vtjocpWoNnbRVtcr+/a9x4MB/CQ5+\nmNatJ1odjssICrJfbOTcW1CQ1ZEpV6IJX9UaycmL2bXrYZo2vZ4OHV6yOhyX4o71eVT104SvaoX0\n9P+xdevNBAVdTLduH2OM/ukqVVH6v0a5vJyc3cTGXo2PTyvCw5fg6elndUhK1Uqa8JVLKyg4TkzM\nMESKiIz8Fh+fZlaHpFStpat0lMuy2fKIi7uG3Nw99Oixgvr1O1sdklK1mvbwlUsSsbFt2x2kpf1C\nt24f0bDhpVaH5NK0lIFyhPbwlUvas2cKR4/Op337aTRvfpPV4bg8LWWgHKE9fOVyDh58j337XqBl\ny3to0+bvVoejlNvQhK9cyvHj37Fjx700bjyUTp3ecvlSx0rVJk5L+MaYp40xB4wxW07ehjmrLeUe\nMjP/ID7+BgICIggL+wwPDx1xVKo6Oft/1HQRednJbSg3kJubREzMcLy8GhIR8Q1eXjrbqFR10y6U\nslxhYTqxscMpKkonOnoN9eq1tjokpdySs8fw7zfGxBhjZhljGjm5LVUL2WwFxMffSHZ2At27f05A\nQKTVISnltqqU8I0xPxpj4kq4jQLeBtoDUcAh4JVS9nGPMWajMWZjcnJyVcJRtYyI8Oef93HixHd0\n7vwujRtfYXVISrm1Kg3piMgQR15njHkP+KaUfcwEZgL06tVLqhKPql327ZvGoUPv067dk7RseafV\n4Sjl9py5SqflGXevBeKc1ZaqfY4cmceePU/QvPkthIT8y+pwlKoTnDlp+x9jTBQgQCLwNye2pWqR\n1NTVbNt2Bw0aDKBr1w90rb1SNcRpCV9Exjlr36r2ysraRlzcNfj5tSc8fDEeHvWsDkmpOkPPtAVY\ntw5efNH+UzlNfv5RYmOHYYw3ERHL8PbWhVtK1SRdh79uHVx+OeTng48PrFgBfftaHZXbKSrKJjb2\navLzDxMV9TN+fqFWh6RUnaM9/FWr7Mm+qMj+c9UqqyNyOyJFbN16KxkZGwgLm09Q0EVWh6RUnaQJ\nf+BAe8/e09P+c+BAqyNyO7t2PcKxY4vp2PE1mjYdZXU4StVZOqTTt699GGfVKnuy1+GcapWU9F+S\nkl4jOPhBgoMfsDocpeo0TfhgT/Ka6KvdsWNfsXPngzRtei0dOmgNPaWspkM6yinS0zeQkDCWwMCL\n6NbtE4zxtDokpeo8Tfiq2uXk7CE2dgQ+PhcQEfE1np71rQ5JKYUmfFXNCgpOEBs7DJECIiKW4ePT\n3OqQlFIn6Ri+qjY2Wx5xcdeSk7ObHj1+wN+/q9UhKaXOoAlfVQsRYdu2u0hL+5lu3ebSsOFlVoek\nlDqHDumoapGYOJWjR+cSGvo8LVrcbHU4SqkSaMJXVXbo0Cz27n2Oli3vpm3bx60ORylVCk34qkqO\nH/+e7dvvoVGjq+jUaYaWOlbKhWnCV5WWmRlDfPwN+Pt3p3v3z/Dw8LY6JKVUGTThq0rJyztAbOxw\nPD2DiIhYipdXkNUhKaXKoat0VIUVFmYQEzOcwsJUoqPX4OsbbHVISikHaMJXFWKzFZKQcCNZWXFE\nRi4lIKCH1SEppRykCV85TET488+JHD++nM6dZ9K48VVWh6SUqgAdw1cO27//Pxw6NJO2bR+nVau/\nWh2OUqqCNOErhxw58im7dz9G8+ZjCQ19zupwlFKVoAlflSs1dQ3bto2nQYP+dO36Icbon41StZH+\nz1Vlys7eTlzcKHx9QwkP/xIPj3pWh6SUqiRN+KpU+fnJxMQMwxhPIiOX4e3d2OqQlFJVoKt0VImK\ninKIixtJfv4hoqJW4ufX3uqQlFJVpAlfnUfExtatt5Ke/j+6d/+coKCLrQ5JKVUNNOGr8+za9SjH\njn1Bhw7TadbsWqvDUUpVEx3DV2dJSnqTpKRXad36Adq0edDqcJRS1UgTvip27NjX7Nw5mSZNRtGx\n46tWh6OUqmaa8BUA6ekbSUgYQ2BgT8LC5mKMp9UhKaWqmSZ8RU5OIrGxI/DxaU5ExNd4evpbHZJS\nygl00raOKyhIJTZ2GCJ5RESsxMenhdUhKaWcRHv4dZyHh/fJK1Ytxt+/m9XhKKWcSHv4dZynpz/d\nuy+0OgylVA3QHr5SStURVUr4xpjRxph4Y4zNGNPrnOceN8bsNMZsN8bolTKUUspiVR3SiQOuA949\n80FjTBgwBugOtAJ+NMZ0FpGiKranlFKqkqrUwxeRrSKyvYSnRgGfikieiOwBdgK9q9KWUkqpqnHW\nGH5rYP8Z95NOPnYeY8w9xpiNxpiNycnJTgpHKaVUuUM6xpgfgQtKeOqfIvJVVQMQkZnATIBevXpJ\nVfenlFKqZOUmfBEZUon9HgDanHE/+ORjSimlLOKsIZ0lwBhjTD1jTCjQCVjvpLaUUko5wIhUfhTF\nGHMt8AbQDEgFtojIVSef+ydwJ1AIPCgi3zqwv2Rgb6UDco6mwDGrg6gBdeE49RjdR104zoocYzsR\naVbei6qU8OsCY8xGEelV/itrt7pwnHqM7qMuHKczjlHPtFVKqTpCE75SStURmvDLN9PqAGpIXThO\nPUb3UReOs9qPUcfwlVKqjtAevlJK1RGa8MtgjBl6strnTmPMY1bHU92MMW2MMSuNMQknq55Otjom\nZzHGeBpjfjfGfGN1LM5ijGlojFlkjNlmjNlqjOlrdUzVzRjz0Mm/1ThjzHxjjK/VMVUHY8wsY8xR\nY0zcGY81Nsb8YIz58+TPRlVtRxN+KYz9Kt5vAX8BwoCxJ6uAupNC4P9EJAzoA0x0w2M8ZTKw1eog\nnOx1YLmIdAV64GbHa4xpDTwA9BKRcMATe1VedzAbGHrOY48BK0SkE7Di5P0q0YRfut7AThHZLSL5\nwKfYq4C6DRE5JCKbT/6egT1BlFjkrjYzxgQDw4H3rY7FWYwxDYDLgA8ARCRfRFKtjcopvAA/Y4wX\nUB84aHE81UJEVgPHz3l4FDDn5O9zgGuq2o4m/NI5XPHTHRhjQoBo4H/WRuIUrwF/B2xWB+JEoUAy\n8OHJoav3jTH+VgdVnUTkAPAysA84BKSJyPfWRuVULUTk0MnfDwMtqrpDTfgKY0wA8Dn2EhjpVsdT\nnYwxI4CjIrLJ6liczAu4EHhbRKKBLKphCMCVnBzDHoX9w60V4G+MudXaqGqG2JdTVnlJpSb80tWJ\nip/GGG/syX6uiHxhdTxO0A8YaYxJxD4sN9gY84m1ITlFEpAkIqe+oS3C/gHgToYAe0QkWUQKgC+A\nSyyOyZmOGGNaApz8ebSqO9SEX7oNQCdjTKgxxgf75NASi2OqVsYYg33Md6uIvGp1PM4gIo+LSLCI\nhGD/N/xJRNyuVygih4H9xpguJx+6HEiwMCRn2Af0McbUP/m3ezluNjF9jiXA+JO/jweqfP2Rql7T\n1m2JSKExZhLwHfbVALNEJN7isKpbP2AcEGuM2XLysSdEZJmFManKux+Ye7KDshu4w+J4qpWI/M8Y\ns9E0giMAAABlSURBVAjYjH2F2e+4yRm3xpj5wECgqTEmCXgKmAZ8Zoy5C3sV4Rur3I6eaauUUnWD\nDukopVQdoQlfKaXqCE34SilVR2jCV0qpOkITvlJK1RGa8JVSqo7QhK+UUnWEJnyllKoj/h9UG6I7\nH99NPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1aecb717940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualizeData(dataArr, labelArr, b, alphas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.3.引入核函数\n",
    "### 6.3.1.在线性SVM的基础上进行改动"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space\n",
    "    m,n = shape(X)\n",
    "    K = mat(zeros((m,1)))\n",
    "    if kTup[0]=='lin': K = X * A.T   #linear kernel\n",
    "    elif kTup[0]=='rbf':\n",
    "        for j in range(m):\n",
    "            deltaRow = X[j,:] - A\n",
    "            K[j] = np.dot(deltaRow, deltaRow.T)\n",
    "        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab\n",
    "    else: raise NameError('Houston We Have a Problem -- \\\n",
    "    That Kernel is not recognized')\n",
    "    return K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "class optStruct:\n",
    "    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters \n",
    "        self.X = dataMatIn\n",
    "        self.labelMat = classLabels\n",
    "        self.C = C\n",
    "        self.tol = toler\n",
    "        self.m = shape(dataMatIn)[0]\n",
    "        self.alphas = mat(zeros((self.m,1)))\n",
    "        self.b = 0\n",
    "        self.eCache = mat(zeros((self.m,2))) #first column is valid flag\n",
    "        self.K = mat(zeros((self.m,self.m)))\n",
    "        for i in range(self.m):\n",
    "            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calcEk(oS, k):\n",
    "    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)\n",
    "    Ek = fXk - float(oS.labelMat[k])\n",
    "    return Ek"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej\n",
    "    maxK = -1; maxDeltaE = 0; Ej = 0\n",
    "    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E\n",
    "    validEcacheList = nonzero(oS.eCache[:,0].A)[0]\n",
    "    if (len(validEcacheList)) > 1:\n",
    "        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E\n",
    "            if k == i: continue #don't calc for i, waste of time\n",
    "            Ek = calcEk(oS, k)\n",
    "            deltaE = abs(Ei - Ek)\n",
    "            if (deltaE > maxDeltaE):\n",
    "                maxK = k; maxDeltaE = deltaE; Ej = Ek\n",
    "        return maxK, Ej\n",
    "    else:   #in this case (first time around) we don't have any valid eCache values\n",
    "        j = selectJrand(i, oS.m)\n",
    "        Ej = calcEk(oS, j)\n",
    "    return j, Ej"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def updateEk(oS, k):#after any alpha has changed update the new value in the cache\n",
    "    Ek = calcEk(oS, k)\n",
    "    oS.eCache[k] = [1,Ek]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "def innerL(i, oS):\n",
    "    Ei = calcEk(oS, i)\n",
    "    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):\n",
    "        j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand\n",
    "        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();\n",
    "        if (oS.labelMat[i] != oS.labelMat[j]):\n",
    "            L = max(0, oS.alphas[j] - oS.alphas[i])\n",
    "            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])\n",
    "        else:\n",
    "            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)\n",
    "            H = min(oS.C, oS.alphas[j] + oS.alphas[i])\n",
    "        if L==H: return 0\n",
    "        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel\n",
    "        if eta >= 0: return 0\n",
    "        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta\n",
    "        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)\n",
    "        updateEk(oS, j) #added this for the Ecache\n",
    "        if (abs(oS.alphas[j] - alphaJold) < 0.00001): return 0\n",
    "        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j\n",
    "        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction\n",
    "        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]\n",
    "        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]\n",
    "        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1\n",
    "        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2\n",
    "        else: oS.b = (b1 + b2)/2.0\n",
    "        return 1\n",
    "    else: return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO\n",
    "    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)\n",
    "    iter = 0\n",
    "    entireSet = True; alphaPairsChanged = 0\n",
    "    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):\n",
    "        alphaPairsChanged = 0\n",
    "        if entireSet:   #go over all\n",
    "            for i in range(oS.m):        \n",
    "                alphaPairsChanged += innerL(i,oS)\n",
    "                #print(\"fullSet, iter: %d i:%d, pairs changed %d\" % (iter,i,alphaPairsChanged))\n",
    "            iter += 1\n",
    "        else:#go over non-bound (railed) alphas\n",
    "            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]\n",
    "            for i in nonBoundIs:\n",
    "                alphaPairsChanged += innerL(i,oS)\n",
    "                #print( \"non-bound, iter: %d i:%d, pairs changed %d\" % (iter,i,alphaPairsChanged))\n",
    "            iter += 1\n",
    "        if entireSet: entireSet = False #toggle entire set loop\n",
    "        elif (alphaPairsChanged == 0): entireSet = True  \n",
    "        #print( \"iteration number: %d\" % iter)\n",
    "    return oS.b,oS.alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calcWs(alphas,dataArr,classLabels):\n",
    "    X = mat(dataArr); labelMat = mat(classLabels).transpose()\n",
    "    m,n = shape(X)\n",
    "    w = zeros((n,1))\n",
    "    for i in range(m):\n",
    "        w += multiply(alphas[i]*labelMat[i],X[i,:].T)\n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3.2.在测试中使用核函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def testRbf(k1=1.3):\n",
    "    dataArr,labelArr = loadDataSet('testSetRBF.txt')\n",
    "    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important\n",
    "    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()\n",
    "    svInd=nonzero(alphas.A>0)[0]\n",
    "    sVs=datMat[svInd] #get matrix of only support vectors\n",
    "    labelSV = labelMat[svInd];\n",
    "    print( \"there are %d Support Vectors\" % shape(sVs)[0])\n",
    "    m,n = shape(datMat)\n",
    "    errorCount = 0\n",
    "    for i in range(m):\n",
    "        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))\n",
    "        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b\n",
    "        if sign(predict)!=sign(labelArr[i]): errorCount += 1\n",
    "    print( \"the training error rate is: %f\" % (float(errorCount)/m))\n",
    "    dataArr,labelArr = loadDataSet('testSetRBF2.txt')\n",
    "    errorCount = 0\n",
    "    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()\n",
    "    m,n = shape(datMat)\n",
    "    for i in range(m):\n",
    "        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))\n",
    "        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b\n",
    "        if sign(predict)!=sign(labelArr[i]): errorCount += 1    \n",
    "    print( \"the test error rate is: %f\" % (float(errorCount)/m) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "there are 30 Support Vectors\n",
      "the training error rate is: 0.130000\n",
      "the test error rate is: 0.150000\n"
     ]
    }
   ],
   "source": [
    "testRbf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def img2vector(filename):\n",
    "    returnVect = np.zeros((1,1024))\n",
    "    fr = open(filename)\n",
    "    for i in range(32):\n",
    "        lineStr = fr.readline()\n",
    "        for j in range(32):\n",
    "            returnVect[0,32*i+j] = np.int(lineStr[j])\n",
    "    return returnVect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def loadImages(dirName):\n",
    "    from os import listdir\n",
    "    hwLabels = []\n",
    "    trainingFileList = listdir(dirName)           #load the training set\n",
    "    m = len(trainingFileList)\n",
    "    trainingMat = np.zeros((m,1024))\n",
    "    for i in range(m):\n",
    "        fileNameStr = trainingFileList[i]\n",
    "        fileStr = fileNameStr.split('.')[0]     #take off .txt\n",
    "        classNumStr = np.int(fileStr.split('_')[0])\n",
    "        if classNumStr == 9: hwLabels.append(-1)\n",
    "        else: hwLabels.append(1)\n",
    "        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))\n",
    "    return trainingMat, hwLabels "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def testDigits(kTup=('rbf', 10)):\n",
    "    dataArr,labelArr = loadImages('trainingDigits')\n",
    "    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)\n",
    "    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()\n",
    "    svInd=nonzero(alphas.A>0)[0]\n",
    "    sVs=datMat[svInd] \n",
    "    labelSV = labelMat[svInd];\n",
    "    print (\"there are %d Support Vectors\" % shape(sVs)[0])\n",
    "    m,n = shape(datMat)\n",
    "    errorCount = 0\n",
    "    for i in range(m):\n",
    "        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)\n",
    "        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b\n",
    "        if sign(predict)!=sign(labelArr[i]): errorCount += 1\n",
    "    print (\"the training error rate is: %f\" % (float(errorCount)/m))\n",
    "    dataArr,labelArr = loadImages('testDigits')\n",
    "    errorCount = 0\n",
    "    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()\n",
    "    m,n = shape(datMat)\n",
    "    for i in range(m):\n",
    "        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)\n",
    "        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b\n",
    "        if sign(predict)!=sign(labelArr[i]): errorCount += 1    \n",
    "    print (\"the test error rate is: %f\" % (float(errorCount)/m) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "there are 61 Support Vectors\n",
      "the training error rate is: 0.000000\n",
      "the test error rate is: 0.010753\n"
     ]
    }
   ],
   "source": [
    "testDigits(('rbf', 20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
